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CHAPTER 1 


INTRODUCTION 

This contract was concerned with die application of recently 
developed function -space Davidon -type techniques to the shuttle accent 
trajectory optimization problem, and with an investigation of the 
recently developed PRAXIS algorithm for parameter optimisation. 

The PRAXIS algorithm has been programmed into the NASA- JSC 
PEACE parameter optimisation program, while the function -space 
algorithms are contained in a separate single program. 

A L the outset of this analysis the major deficiency of the function- 
space algorithms was their potential storage problems. Since most 
previous analyses of the methods were with relatively low-dimension 
problems, no storage problems were encountered. However, in 
shuttle trajectory optimisation, storage is a problem and this problem 
was handled effectively. This point will be discussed further in Chapter 4. 

la Chapter 2 the function -space algorithms are presented and 
discussed. The theory is presented in such a way that both parameter 
and function controls are handled naturally. Aerospace problems 
usually contain both types of controls. 

In Chapter 3 the shuttle ascent model is presented along with the 
development of the particular optimization equations, hi Chapter 4 
the operation of the algorithm and typical simulations are presented. 
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fa Chapter 5 varia le final-‘ime problem considerations are 

studied since some investigators have found the fun ctioi. -space 

accelerated- gradient methods to behave poorly on variable final-time 

problems. Simulations and heuristic reasoning indicate that the initial 

( 0 ) 

choice of t^ (in the iteration scheme), say t^ , is critical, and that 

( 0 ) - 

a ch o i ce t^ < t^ (optimal t^) appears to improve the convergence 
rate considerably. 

fa Chapter 6 a modification of Powell' s algorithm developed 
2$ 

by Brent , is presented and discussed. The algorithm is known as 
PRAXIS, and it is a parameter optimization scheme which does not 
require gradient information. A flow diagram of the algorithm is 
presented in Appendix D, and a listing, is available with the NASA- TSC 
PEACE program. Finally, Chapter 7 presezt? the conclusions and 
recommendations for further study. 



CHAPTER 2 


THE ALGORITHMS 

b the past few years numerous quasi-Newton type algorithms for 
the solution of parameter optimisation problems have beta extended from 
Euclidean spaces to infinite dimensional real Hilbert spaces. Just as in 
Euclidean space, the primary advantage in Hilbert space is the accelerated 
rate of convergence due to file building of second order information while 
requiring only function and gradient evaluations. Except for the conjugate 
gradient and gradient methods, existing function space methods cannot 
handle directly control variable inequality constraints. Thus applications 
to optimal control problems have primarily dealt with the classical Bolsa 
problem. Since most realistic problems cor tain control variable inequality 
constraints it is desirable to be able to handle them directly in a computation 
scheme, hi attempting to solve such problems a new function space 
algorithm has been generated and two existing quasi-Newton type algorithms 
have been modified to allow them to handie directly the bounded control 
problem. The modification of the algorithms was strongly influenced by the 
work c £ Pagurek and Woods ide^ on extending the conjugate gradient method to 
t&elude bounded controls. The methods modified include Davidon and 

3 

Broyden type algorithms. 


3 
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2.1 The Algorithms 

Is this section the various algorithms will be formally stated for 

the problem of minimising a real functional J(u), where u may be either finite - 

or infinite -dimensional. With u finite -dimensional, the formulas are 

applicable to the standard unconstrained parameter optimisation problem. 

In tbe next section, 'Jte appropriate modifications for application to optimal 

control problems with bounded control variables will be presented. 

In die listing below, each algorithm requires the specification of 

a starting vector, u . In addition, the Davtdon and Broyden 

o 

algorithms require the specification of a positive -definite, self-adjoint 
linear operator, H q . Also, <a , b> and a><b will be used to denote the inner 
and outer (dyadic) products, respectively, on the given Hilbert space. 

T 

Note that if the space is n-dimensional Euclidean space, then <a, b> = a b 
T 

and a> <b - a b , where a, b are n-dimensional vectors. The inner and 
outer products for the optimal control problem will be defined in the next 
section. 

Let g(u) denote the gradient of J, and define the update formula by 

u = u. + **. d., (2-1) 

i+I i li 

where d^ * search lirection vector and a. 5 scalar parameter defined 
by a one -dimensional search technique which minimizes J with respect to «. 
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I. Gradient Algorithm (G) 

a. Calculate the search direction d. - -g(u.). 

i 1 

b. Use Eq. (2.1) to calculate u. + ^ and return to step a. 

II. Conjugate Gradient Algorithm 1 (Cd)^ 


a. Calculate the search direction* 


( 2 . 2 ) 


where 


P 


V 

‘- 1 <* 1 . 1 - h-f 


(2.3) 


b. Use Eq. (2.1) to calculate u.^ and return to step a. 
m. Davidon Algorithm (PAY) 2, 7 

a. Calculate the search direction 

d. = -H. g 'u. ) (2.4) 

i i ° i 

b. Use Eq. (2. 1) to calculate u.^ 

c. Calculate 


s. = u. . -u. 
i t+1 i 


n * « ,u i Tl > - g(,1 i> 

d. Update H according to the following formula; 


s.xs 

H. .. = H. + — — 

i+l i <s., y.> 


H.y.xH y. 

11 4 1 

<y., H.y.> 

' l ii 


(2.5) 

( 2 . 6 ) 


(2.7) 


e. Return to step a. 


* On the first iterate ( i = 0), define = -g(u^). 
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The same as DAV except for step d, where H is updated according 




s.xH.y. 

i XI 


H.y.x s. 

li x (2.8) 


< s., y.> < s., y.> < s. t y.> 


i 



j ; :x*' r *^ s * o*~< lo v.ontrol Frg-'j^nj; 


As anted previously, in n -dimensional space the algorithms are 

nsed to minimize a scalar valued function .T (u), where u is an n dimensional 

— T 

vector , tiie inner product is < s. y> — s y, the dyadic operator is 

s X y ~ sy X , and the H operator is an n x n matrix of scalars. 

Implementation of the algorithms on this type of problem is well documen*ed 

in the literature. All of the algorithms described, with the exception of tb 

(BRD) algorithm, have als_ ’>een gt cralized to optimal control problems 

where g is the gradient of a functional. The primary difficulty iu implement: 

the quasi -Newton type algorithms on optimal control problems lies in 

reneesen.Hng the infinite -dimensional H-operator. . 

-/ T 

In L., mace the inner product is < s, y> ~ f s ~ v dt and the 

— 8 9 fc r> 

dyadic operator is (sxy) a ~ <y, u> s ’ . However, there simply is 
no convenient way to represent H. One way to overcome this difficulty is 
preserved in Reference 4 by Las don, where it is observed that only H.p^ 

•'rot H. itself) is needed to compute d.. This is also true for the Broyden 
a'ccrithm. To implement the Broyden algorithm, where g is the gradient 
'ho frsr. '".oral, and u, s, and v are time functions, we proceed as fellev 
i: is taken to be any positive -definite, self-adjoirt operator. 

ii> F.m~e=,s H. in Fq. f2.-' as a srro br ’?• II . C'-'err*'' or. »he 

1 A .-N 


rr expression for It. with <*. f 


c < rch 
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i-1 


d i'- H o*i-jS 


L , <Y i- Vi > ) < V«« > . *%Z± il . 

\ < v y ) < v y j > s < vv 1 


*V «!» 

<s j- y 


H. yj 


(2.9) 


Equation ( 2.9 ) requires tike computation of inner products of tike 

functions H.y., s., and y. f and operating with H (H = I being the simplest 
ill l o o 

choice). The functions (s^, .... s. ^ ) are available from past iterations. 
To compute the functions H.y.» we need only replace -g^ by in Eq^2. 9 )» 
i. e. , H. operating on y. instead of -g^. Then, for the case i-1: 

H y = H y + 2 j ( j + Hilfilil) HlJirll , 

w i-l y i-l o y i-l \y + <s.. y> ) < b : , y.> a j 




< s. 


r i-r 


<*•» y > 
J y j 


J 


< W 


H.y. 
J J 


( 2 . 10 ) 


Thus Hj j 7^ j can be computed in a way requiring only inner products and 
operation with H = I, as was the case for -H.g.. Note that 2i + 4 time 

O 11 

functions must be stored after the i iteration in order to compute the i+1 iterate, 
‘ * i+1 functions 

i functions (2.11) 


(»o’-*" ®i> 


<V, H i-i Vi> 


Vr Vi 


3 functions 
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We shall now define the basic optimal control problem, and then discuss 
the problems of implementing the quasi-Newton algorithms. The 
interpretation of the above formulas and operations is more motivating 
in an optimal co: .rol setting. 

The optimal control problem of interest is a Bolza problem with 
control constraints as follows: 


Minimise:. 

J(u) = +(x { ) + / L (t,x,u) dt 

t 

o 

(2.12) 

Subject to: 

x = f(t, X, y), x(t ) = X (x = n-vec jt ) 
o o 

(2.13) 


Jllj | <_ (l = If • • * >Btl) 

(2.14) 


t Q , tj specified 

Terminal conditions are included in the $ (x^) - term and statevariable 
inequality Constraints are included in the L( t, x, -u) by the method of 
penalty functions. 

A motivating way of viewing the -quasi-Newton methods is as a 
class of algorithms between the first-order*® and second-order *®’ *** ^ 

optimal control gradient methods . The goal of a quasi-Newton algorithm 
is to build information about the second -variation operator without computing 
it explicitly, i.e., based upon gradient information only. Since only 
gradients and function evaluations are required for the quasi -Newton 
methods, we shall first outline the gradient method for optimal control 
problems, and then discuss the modifications for a quasi -Newton method. 
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In all of the algorithms, the following equations are required: 


H = i. J \ fit. t . u) 

X=-*I Xft)-^ 

a* ' ax. 


. . an 
« w * ^ 


(2.15) 

(2.16) 
(217) 


The function H above is the Hamiltonian, which is not to be confused with 
the oper a tor H. of the algorithms, and g(u) = 8H/8u is the function space 
gradient. The usual implementation of the standard gradient method is 
shown in Figure 1. where 

"in = V (2-i8) 

X 

Note that the subscripts indicate Hie iterate number for the respective 
vectors; this allows less cumbersome writing of the quasi-Newton formulas. 

The optimal control f r- the problem defined by Eqs. (2.14-16 ) will, 
in general, consist of a sequence of interior ( |u_ j < c.) fad bounded 
( ?u. { = c. ) control component intervals. On each subarc Hie following 


renditions most be satisfied: 


u - c -- 1 — 
l 1 

m 

X 

< 0 

(2.1?) 

•c. < u, < c. " 
111 


&. =o 

Pu. 

i 

(2.201 

M C — 

k- 

m 

e-x. 

^ 0 

(? "’M 


i 
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We shell now discuss how bounded control variables ere treated directly 
in die standard gradient method since die same basic idea is employed in 
die quasi -Newton methods. 

As new controls are generated by varying a in Eq. (2.18), they may 
violate |u. i < c.. On these intervals a is truncated such that if u. > c., 

1 l'“ i it 

u. is set equal to c. andif u_< -c. , u. is set equal to -c_. After truncation 

t i til t 

die cost associated with die gives a is calculated, hi this way die saturation 
region may change from iter»tc co iterate and costs are only computed for 
realixeable controls. 

The i mp l em entation of die quasi-Newton type algorithm on unbounded 
control problems it shown in Figures 2 and 3. As in Figure 1 die subscripts 
indicate the iterate number, and Eq. (2.7) implies die Davidon algorithm 
and Eqs. (2.9) and (2.10) imply the Broyden algorithm. 

Note that as the iteration proceeds, the number of functions stored 
increases. The computation time per iteration will also increase because 
of more inner product evaluations in die updating formulas for Hy and d. 

To overcome this difficulty die algorithm is restarted with m pure gradient 
step when i = q, where q is some predetermined interger. Pierson and 
Rajtora** demonstrated that the restart feature sometimes speeds 
convergence in addition to being a practical necessity. For certain problems, 
they found that for q small, say 3 or 1, that the convergence rate was 


enhanced. 



1 * 






u. 


— Tota 1 sf ."age g Tf u.. Sf> , y pt Ho p 

s ! 

O'.v of the Function S:>u*o Quasi- 7 "c\vton /.Iporithms for 

H - I .-pi i - 0, 1 

o 


Figure 2. 



i*n i"2 


»>(t) 


*,<t) 


! ** X 2»f> = ♦- 


iCompute X 2 (t), g 2 (t) = H ft) 


Store y|(t) = g 2 (t) - gj(t) 


V 

. Calculate and store Hjy. = Iy^ + £ [ B%'. 2»7 or 2. 8 ] 

j=0 

1 

Calculate ^ = -g 2 (t) - 2^ [ Eq. 7 or 8 ] 


1-D search => ». 
— U 3 = tt 2 + * 2 d 2 


Store s 2 (t) = u 3 (t) - u^t) = « 2 d 2 

Total storage g 2 , Hy s Q , j v 


i s** W = ♦, 


p-npote SiW. g„tt) * H„(t) 


Stor ® Vl W = 8 ^ <,, - «n-l W 


n-6 

— Calculate and store H .y . = Iy , +11 [Eq. 2.7 or 2.8 ] 

n-ln-1 n-1 j=o 

n-1 

''ilculate d^ = ~g n (*) - [ Eq. 2.7 or 2. 8 ] 

1-D search sOft 


- u . = u + <* d 
n+1 n n n 


Star. s n (t) = u M1 W - u a (t) = « n d n 


Total storage g^ « . . T B .,. Vo 


H .y . 
n-ln-1 


Figure 3. How of the Function Space Quasi-Newton Algorithms for 

H — I and i — 2 p . . . » n g« « . 
o 




2. 3 Bounded Controls 
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To apply tiie quasi-Newton type algorithms to the bounded control 
problem a modification to the updating formula is required. In the interior 
portion of the control we wish to build second order information while 
second order information on the bounded portion of the control is of little 
use. Thus, tiie quasi-Newton formulas should concentrate on the interior 
controls, and a standard gradient formula can be usei on 'he bounded 
portion of the control. 

As with the gradient algorithm, as new controls are generated 
they are truncated before calculating the associated cost. A 
saturation function w.(t) identical to Pagurek and Wood side's* is defined. 
This saturation function is set equal to sero when the control is on the 
boundary and is set equal to unity on the interior. The saturation function 
is then used in the following way to compensate for our lack of freedom 
in choosing the control on the boundary. Instead of using q. y» and Hy in 
the formulas for calculating the search direction and updating Hy. we use 
wg. wy, and wily. We know that g = 0 on the interior portion of the 
optimal control and this is where we wish to build second order information. 
On the region of saturation g ^ 0 (in general) and the y's (or Ag's) should 
not contribute to the inner products in the updating formulas. It is not 
necessary to apply the saturation function to r because on the saturation 
region, s = Au will already be sero. 
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The quasi -Newton algorithm for bounded control problems differs 
from the algorithm for unbounded control problems in the following ways, 

i) As the 1-D search seeks die best a die associated controls are truicat*- 
before the associated cost is calculated. 

ii) A saturation function w.(t) is generated after each iteration. 

iii) wg, wy, and wHy are used in the updating formulas and in the 
c-)culation of the search direction. 
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2.4 Function and Parameter Controls 

Some optimisation problems are most naturally formulated 

using a combination of controls ui R n and spaces. The shuttle 

ascent problem developed in the following chapter is such a problem. 

While the approach for problems whose control space lies in either 
n n 

R or L is well develr ped, the theory is incomplete when a mixture 

€m 

of controls exists. For die general Balsa problem Eq. 2.12-14, the 
control is an. M vector of functions. 


^(t) 

l 

I 


' • • f V V 


and the first variation after appropriate adjoint function definitions is, 

tf T 

H 6u dt (2.22) 

„ u 

'"O 


« J = 


-l 


The gradient of J at the element u, denoted by g, is defined by the inner 
product relation, 

5J=|- J[S(t)+«h(t)] | f =Q • = <g. h>« =<g,6u> (2.23) 

On I® f t Q . t £ ] <v, s> =/ t V T s dt (2 . 24) 

■ o 

Comparing Eqs . 2.22 and 2. 23 we have, 

g(t) = H u (2.25) 

This is the usual function space gradient, ft can be shown that 
in L.® [ t Q , t £ ] the linear quadratic problem ( L Q P), 
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min J = 7 X [ * T P <*) x + u T R (t) u] dt 

b V 


SUBJECT TO x = G (t) x + B (t) a (2.26) 

* ( fc J * * * » t » V « iven 

o o oof 

in equivalent to the minimisation of the unconstrained quadratic functional, 

J = ^- <u, A u > + < u, w > + J (2.27) 

Z o 

where A, w, and J q are appropriately defined and the inner product 
is defined by Eq. 2.24. It can also be shown that A is a strongly 
positive operator if P (t) > 0 and R (t) > Oon [ t , t { ] . fhen, the quasi- 
Newton methods can be applied directly to the quadratic functional 
given in Eq. 2.27 for which convergence can be shown. Therefore an 
iterative solution to the optimal control problem of Eq. 2.26 can be 
generated. 


hi general, the nonquadratic functional is of interest. However, 

if the general Bolza problem can be approximated by a second-order 

expansion in the neighborhood of the minimizing control reasonable 

convergence may occur *near * the minimum. 

Accepting the desirability of the approach above, consider the 

class of optimal control problems whose control space is composed 

of elements in L.^ 11 [ t , t 1 and R n . 

2 1 o f J 

^ (t) 

I 
I 
I 


C 1 

I 

i 

cl 

n 


— . , m ® 

= u * x R 


(2.28) 



Where 


{“l ttK * 

fv — • '„} 


_ in r , 

L Z t V J 


The first variation becomes, 

t t T 

5J = J H 5ttdt+| 6 c dt 

u * t 

o o 

Since c. • constant, then 6 c • dc and Eq. 2.27 becomes 


( 2 . 


6 J = / £ H u T 6 a dt + dc T / t 1 Hf. dt 


(2. 


m 


Proceeding as in [ t^, t^] an inner product must be defined wh 
will imply the gradient. Then to justify, in a sense, the application c: 
the quasi -Newton algorithms to this class of problems a suitable 
optimal control problem similar to the LOP must be developed which 
can be reduced to an equivalent unconstrained quadratic functional. 

The merit of the chosen inner product is measured by its usefulness 
in implementing die quasi -Newton methods . 

A more straightforward approach to this class of problems is 
based on the observation that bounded elements in R n can be chought 
of as constant functions in [ t Q , t^] . Thus the control can be 


partitioned, 


■Y“ 

i 

r 

“p<‘> 


U =: 


U 


P A i 


U 





2C 


where u^ ( i - 1, . . . , p) e L^ [ t©, t^] and u. ( i = p + 1, . . , m) 
are finite constant functions. Thus the admissible control space is a 
subspace S of L™ [ t^, t f ] , 

S = | u | u. ( i = 1, . . , p) c L> 2 [ t Q . t f ] ; ( i = p + 1, . . . , m) 

* finite constant functions (2. 32) 

The goal is to find an* S which minimize s the quadratic 

functional of Eq. 2.27. The advantage of this approach is that all the 

theory which has been developed for L™ [ t^, t^ ] still applies. The 

only change is that u is restricted to lie in S. 

The set S is a linear subspace of L m [ t , t ] . The usefalness 

2 ox 

of the fact that S is a linear sabs pace of L m ft . t 1 lies in the 

2 of 

following property of Hilbert spaces and quasi-Newton algorithms. 

Property : Let M be a linear subspace of a Hilbert space D. 

There exists a mapping P:D <4 M called a 

projection operator which is linear, seliadjoint, 

and indempotent. If of the quasi-Newton 

algorithms is chosen to be P and u « M , then 

o 

u. c m for all i t and 

i 

lin ||H o g k || 2 =0. 
k-* 

that is, the project) rn of fhe gradient ontoM tends 
*o which is - for convergence. 
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it) I dt 


- i ; r *f 

: A i E f •- f \ 

t ' t 

o o 

A. 

r f A t \ d« +f^/A c l dtl|/ tf B c dt) 


TrB.i 

I _ ^ 


r f r A 1 r £ A .* u f 

{ lA t it/ V d ‘l [ r c j 


= <PA, B> 
? 

iii) Idempotent (P = P) 

Ta. if A 


A 


f 


1 

A t J 

f A d 

At 

'■ c 

f 

- 

L 

o 


P = P 


1 r tf a ^ 1 r tf 1 r tf a ^ JT . 

m { A c dt A c dt dr 

o o o 


j- f { A dt 
At - c 

o 


= P -- f - =>P 2 = P 

A 


The properties developed above imply how the first variation of 

La. (2;. 30) should be treated in the quasi-Newton algorithms. First, 

viewing (u. (t), u (t), c, , . . c ) as an element cf L_ m+n t , t,, the 
I m 1 n dot 


"ratlierif is 

r« 


- U 
a r 1 


(2. 34) 
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and an admissible choice for H , ?ay H , is the projection operator (< 

o o 

which implies that the initial search direction is 


~ /V 

H , g = 

C o 


H 


( 0 ) 


a 


t 


f H *°'dt 
- c 


(i. *5) 


However, n* that this is equivalent to assuming 


r 


g = 


H 

a 

'tr'Tov 


i ff 1 i 
ITT / H c * 


f O t 


J 


(2- 36) 


with ( u^t), .... u (t), c^, . . . , c^) « L^jt^, t f |x R n , and 


H g =| 
o o 1 


H 


( 0 ) 


u 


1 t (0) 

t, t f H dt 
f o c 

t 


= H g . 
o o 


(2. 37) 


Furtheiracrc, the choice of definition for the gradient (2 36 i has ►he 
^rne convergence properties as the choice (2. 34) since 
IlSKII-o 


* r 


i » * j » : 

r ,. ■ ' - ' ’ 1 


t t I ’ • n [ 2 . *P ' 

r i 

v. ]\] ho ;o'T c ’ho gradient r\;;;cs f sion i T * Char 


.33), 


I since 



CHAPTER 3 


SPACE SHUTTLE ASCENT MODEL AND OPTIMIZATION 
3.1 Vehicle and Mission Description 

The vehicle and mission considered are taken from Reference 15. 

The goal is to determine a control history for the pressure-fed series 
barn booster/040 c orbiter which will yield maximum payload deliverable 
to a 50 x 100 mn. orbit inclined 28. 5 degrees. The trajectory is constrain "d 
to 650 psf maximum dynamic pressure and 3.0 g maximum accleration. 

For trajectory purposes the mass of the vehicle car be broken down 
into five main subdivisions . 


f MASS DISTRIBUTION i 

| BOOSTER 

ORBITER 




mg j = fuel first stage = 3. 50680 x 10^ lbm. 

m^ = fuel second stage = 1.16415 x 10^ lbm. 
m # = structure first stage = 5.70850 x 10^ lbm. 

5 

m^ ^ - structure second stage = 2.61300 x 10 lbm. 
m^ = payload = quanity to be maximized 

The trajectory is determined by two controls, the mass flow rate 
which implies the magnitude of the thrust and a thrust angle. The mass 
flow rate may vary from zero to ; .laximui' whir’ trust between 

0 and 100%. The overall trajectory can be t. . ^ four "phases". 
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Each of the phases is characterized by the way in which the thrust angle 
is determined and by die coordinate system in which the equations of 
motion are being integrated. 


FIRST STAGE 

■vjj 1 E T>\ ^ 

COORDINATE SYSTEM SPHERICAL ROTATING 

APPENDIX A 

POLAR 

PHASE 1 

PHASE 2 

PHASE 3 

PHASE 4 

VERTICAL 

RISE 

PITCH 

OVER 

GRAVITY 

TURN 

LINEAR 

TANGENT 


The equations of motion for die first stage are integrated in a 


spherical coordinate system which rotates with die earth. This 


coordinate system was chosen because of die ease of representing 
initial conditions and aerodynamic forces. The general equations of 
motion are derived in Appendix A. Assuming die first stage engines are 
perfectly expanded to vacuum pressure we have. 


where 


Thrust = T = I |m I - P A 

sp^ 1 1 atm exit 

I =270.7 sec. 

8 Pl 


A = 700 ft 
exit 


0 < |m | < 3. 01385 x 10 4 lbm/ sec. 


T = 8.15849 x10 lbf 
max 


The first stage burn is further divided into three phases. They are. 


i) Phase 1 - vertical rise for ten seconds 


ii) Phase 2 - pitch over at a constant rate for ten seconds 

iii) Phase 3 - gravity turn i.e. , the thrust is parallel to the velocity. 
This phase terminates when all fuel is exhausted in the first stage. 
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Aerodynamic forces are on the order of 2% of die total forces 
acting on the vehicle after staging and drop oft rapidly. Thus, aerodynamic 
forces are neglected daring the second stage barn. This allows the 
equa tio ns of motion to be integrated in a polar coordinate system. 

- The equations of motion are stated in Section 3. 4. This change 

of coordinate systems results in a new set of state variables. The equations 

relating the state variables before and aftev staging are derived in Appendix 

C. By integrating the equations in a polar coordinate system we have 

reduced the number of state variables, simplified the terminal boundary 

co nd i tio ns, and simplified the adjoint equations . It is assumed that die 

second stage engines are perfectly expanded to vacuum pressure, thus. 

Thrust = T = I [ml, 

® P 2 

where I = 456. 5 sec. 

8P 2 

0 < |m | < 3. 0887 x 10* lbm/sec. 

=> T = 1.40999 xlO 6 lbf. 
max 

During second stage burn toe thrust is orientated according to 
the linear tangent steering law, i.e. , 

tan y = at 4 b, (a, b constants ) ( 3. 1) 

where y is the angle between thrust vector and the local horisontal. 

The above discussion leads to the following overall problem: 

# 

i) Initial conditions - launch from KSC 

ii) Terminal conditions - 50 x 100 nm orbit inclined 28. 5 degrees 
with insertion at perigee. 

iii) Controls and Unknown Parameters 
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1) Init ia l GLOW - a parameter. 

2) Mass flow rate m(t) or, equivalently, thrust a .function of time. 

3) y during ptich over - a parameter. 

4) 4* daring ptich over - a parameter control for out -of -plane 
thrust (explained in Section 3.2). 

5) a and b - parameters used during linear tangent steering 
tan y = a t + b. 

iv) Constraints 

1) Q < 650 psf 

max — r 

2) Acceleration < 3.0 g*s. 

max— * 



J 


Figure 3.1 Coordinate System for First -Stage Computation 
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i) Yc’-tirf.l. Rise: 


X = |Tj o (10 seconds) 


ii) Fitch Over: Consider the trial of vectors e , e_, e. . 

r o f 

The nlane defined bv e« and e is the local horizon il plane. 

^ 9 

The caii vector c points in the easterly direction for® f* 0 or v. 

n 

Alter vertical rise, the vehicle pitches over and at the same 
time the plane of the orbit is determined by thrusting at some 
constant azimuth angle $ . (See Fig. 3.2) 



Figure 3.2 Thrust angles. 


The initial thrust is in the vertical direction, i . e- , y = — - The 
then pitches over with y = constant, which implies 

- Y <t-t Tr >, <t„< t < 20) 

whe,*c t = time of pitch -over initiation. The angle »!» = Constant 
throughout vertical rire, and it is noted that tp will not correspond to the 
final inclination. However, the final inclination will be very strongly 
influenced by ^ and, in fact, 4i will be the primary control which affects 


The vehicle 


the final inclination angle. Thus, 
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It 


I jjSin -y e y 


4- Cos y Sinf + Cos y Cos 


* ej (3.2) 


pT | |^Cos Y « r + Sin y (t-t^ ) ^Sin * t Cos * e 4] 


iii) Gravity Turn: T 1 1 V 


T = |t| — 

!▼! 


( 3 . 3 ) 


T = |f| 7 0+ ^ 7. 

|v| r fv| 8 \v\ 4 


( 3 . 4 ) 


iv) Linear Tangent: tan y = at + b 


Sin y = 

T vl+tan 2 


(3.5) 


(-*/2 < y < */2) 


Cos y = 


*4 fan 2TJ p 


* ♦ 1 


T 



To' K< 


! 

-j 
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3. 3 Aerwiysamic Forces 

B is assumed that the vehicle is aligned with die local wind 
Telocity, thus 


A = Drag = D - V 


|D| P M 2 AC d 

C s C_ (altitude, mach number) I V| = •<] 2 2 2 

D D 1 1 u + v -t w 


A = |A| ^ 

|V« 


— 1 2 2 2 
A = p AC D (u +t^ +w 4 ) 



e 



% 


= -r* p a c nT 


2 2 2 — — — 
u + v + w ue + ve^ w e 


0 


(3. 8) 


In all equations that follow the control notation will be, 
u - mass flow rate magnitude 
c_- GLOW (Gross Liftoff Weight) 


c - pitch-over rate during pitch-cver phase 



n 


Cj - a (linear tangent parameter) 

c^- b (linear tangent parameter) 

Cj- oat -of -plane thrust angle during pitch over. 
3. 4 Equations of Motion 


(b) 


1 3 


(«) * 2 -* 4 /(v R o) 


<«> 


*3 * ( X 4 Z + X 5 Z ) *(*l * R o) ' k/ ( *1 * + ( *1 + R o) ^ S " 


+ 2Q x Sin x t 

5 2 x 6 


(v) 


*3* 


4 ( *1 + R o) **“ *2 ( X l + a .) 


— + ( X. 4- R ^ 0? Sin x_ Cos 

\ VI o) 2 


+ 2 H x_ Cos X- + — 
5 2 x. 


(3.9) 


(w) x 5 = - 


X 3 X 5 


X 4 X 5 


( V-\) ( *i + R o) tar - 


- 2ftt^ Cos x^ 


2x, Q Sin x_ + — - 
3 2 x. 
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(mass) x ^ = -u 

i ) Vertical Rise: Fj = jf | - Q ^ x^ x ? ) C D ( ^ 



F 


3 


= 0 


2 = 



2 


x 


3 


|T| =ul sj) - Ptx^Ag (3.10) 

ii) Picch-Over: F x = }T| Cos C 2 ^ fc-t^ 

-O ( Xj.X3.X4.X5) C D ( V * 3 ,x 4 .x 5 ) x 3 

F 2 = |i|SinC 2 (.-. vr ) SmC 5 

-o( Xj. x 3> X4. C D Xj.Xj.X4, x s) x 4 (3U > 

F 3 = |T| Sin C 2 { t-. vr ) Co. C 5 
' Q x l' x 3’ x 4' x s) C D ( x l' x 3’ x 4' x s) x 5 


I CM 
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^ ^ r 
iv) Linear Tangent 'f'. - |T | i — r — 

1 V r +l 


f 2 - |x| 


Tan T =C t + C 


|T| = u I 


sp 


Equations of Motion 

» 

*1 = x z 

^ = ^r 2 g f v 
2 3 


/ (v R o)- k & < R 0 y * r i h 

*3 ' V ( V R o ) + 1 2 ,5 4 


= -u 


(3.13) 


(3.14) 


3.5 The First Variation 

In order to apply optimization theory the problem must be stated 
In control notation or format. The rc are five parameter - type controls 
and one function- type control. Recall the Parameter Controls: 
ij Cj - GLOW 

ii) C 2 - y 

iii) C 3 - a 

iv) C 4 - b 

v) C 5 -«l» 

Function Control: u (mass flow rate magnitude). The equations of 
motion for the system have already been derived (Appendix A ) and 
may be symbolized by 

x = f (t,x,n) ( 0< t<t g ) 

ir = T( t.x>) 


(»„ < t i t f > 
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where, for convenience, u denotes the vector ( Cj, ... , , u). The 

terminal boundary conditions will be handled by the method of quadratic 
penalty functions. The state variable inequality constraints will be 
handled by integral quadratic penalty functions. The performance index is, 

+ P 2 ( Vf> - *2f ) 2 

+ P 3 ( x 3 ' *3f ) 

‘ s , 2 

+ j ( q-650) U (q-650) dt 


+ P / S (acc-3.05) 2 U (acc-3.05) dt 

5 to 

+ P , f Wc-3.00) 2 U (acc-3.00) dt 

6 t 

s 

+ P ? ( Cos*(t f ) - Cos* f ) 2 


Where 4> - inclination of orbit 


U(t,) = 


(m >e 

(0 r\ < 0 


(3.15) 


q - dynamic pressure 
acc - axial acceleration 


P. - penalty weighting factors 

t - defined by fuel exhaustion 1st stage 
8 


t^ - defined by fuel exhaustion 2nd stage 


( 2 ) 



(4) 


1) Vertical rise 

2) Pitch over ^ £ l°/sec. 

3) Gravity Turn T | | V 


MiriWh., 4 ) Linear Tangent tan y = C t + C 

Figure 3. Phasing of Shuttle Ascent Trajectory Optimization. 3 4 
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The problem reduces to: 

10 - 20 - t - 

Min I (u) = <*> ( x , x g , x.) + f L. ( t,x, u) dt + f L^(t f x, u) dt + ^* S L_(t, x, u) dt 

o 10+ 20+ 


4 LX*. u ) dt 


( 3 . 16 ) 


Subject to: x=f{t, x, u) (t <t<t ) 

O — s 

and / x’= f (t, xT, u) (t < t < t.) 

s — I 


to = 0 sec 


x^ -Xg fixed 
x^ free 


10 sec i 20 sec 


t free i t free 
s I s 


| t^ free 


x -x free i x(t )?§ (x(c-))| xT-xl free 
13 , s s i I 3 

»I» (x, ) =0 ' i ib (x. ) =0 

s V t» ' i f 4 ' 

„ ' i 

\ i 

.Mass of Fuel > Trans- J Mass of fuel 

1 st stage i formation | 2 nd stage 

depleted de- j Equations [ depleted defines 

fines t i App. C ! t, 

s , - ^ ! f 


Define 


H = L*h\ f(on each subarc) 


In- first variation is 


6 J = 


. i , T T ^ 

9 dx t <J> dx 4" dx r 

xc o xs s xf f 

+ ( [l^S x + H T &u- \ T 6 i] dt 
0 x 


+ f [H^ 5 x + K^ 6 u - V^ 6 x] dt 


+ [H (t-) - X (t-) x U \V 
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b f 6 5 x +- t-i^S u - fix] dt 

20 + ^ 


u 


j. /\-T '>w + — 

-[fity- 1 (t+i x n 3 )] ^ dt s 
+ [h (t { ) -t T iljV-Sf. «t f Sf «t f 


t. 

I 


+ / [k 1 sx -i-h *u :^wr.ldi 

•i x u 


A.I. 


(3-17X 


T_ . 


Then, integrating -X fix by parts 

, 'T 1 

'ta 


X/ -X A 6i: dt = X* (t 3 ) fix (t a ) -\ T (y fix (y + X a fc \ T 6x dt 


and substituting into Eq. (3,17) 

5 J = <f> ^dx + <t> ^dx + dxl 

xo o xs XI f 


10 


T t * vf *i* <p 

+ X. (o) fi X ( 0 ) -X CIO) 6 X (lol + r [(H :- \ ) A fix i H , ou] dt 

~o x ' 

+ X T (10) fix (10) -X T (20") 6 x(20> f [( H + X) r fix + H 1 fiu] dt 

io x u 


+ [H(t : ) -X (t-) x (t-)]dt 


+ X. (20) 6x (2(f) (t- ) 6x(t-)+ f fHx +*1) 6x + H 8ul dt 

S 8 20 U 


( 3 . 18 ) 


- |h (t+) - / X T (t+) x" <t+>] dt 
l* S S SJ S 

+ |li(t f ) 3: T (t f )‘Sr(t f )J dt { 

+")? (t+) fijT(t+) (t ) 6x"(t > + f |(H ^T) T fx*t H fiuldt 
s s i t • l U J 
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At t *nd need to substitute in the relationship between the variation 

it* x end the t*iff-*rtiai in » 

rV - •*■ >- r -f-f 

*r - ^ - v <i; ( 3 . VM 

Af^er bnhsritntinn end reduction, *he follovinp cxnrf ss'fp is ob^arro.d 

dx„ 

f\ -X It- 'I' d- «- T T |til dv' 
t X .s S 

L - - l f * 

+ [h ft—) -H(£*-i| dt 
L s sj s 

HfV* f 

f r ' t i 1 

* f I - v r.x TT *f -,| dr 

V. *- : * u J 

ZC r t - 

+■ f I ,1X + X» SstH 1 S n| d* 

lOvl x n J 


/ c »“ * t — 1 

^ f "jf^ T t Sv i 1 

% *> »j- ■ % 

«# * « 

-*- f TtyT - v\ f xj n I d* 
I -• ~ n J 


^ ,ff C!Tcr*uds fix ^ ]>r **>r> o r ' ?> 


1 n t -» -r* q f n" — r r» * • ~ ' 


T- -]. 


dx 


1 


•-•■o Eq. < ‘ 
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[♦ -X(t-)] T dxU^McS' = 0 

a s s s s 

s 



[\ T - ^‘a-S + * T <‘J> fH. ] ^ =o 

f s 

Then, since dx g is arbitrary 

X(« s -,= X(«J) |f I 

* 8 


(3. 21) 


(3.22) 


Finally, 


61 ‘ [\o * X ‘ ‘ Vi» ^ ‘v] dm o 

10- — 20- t- t 

+ f H T 6u dt + f H T 6u dt + f S H T 6udt + f f hF 

0 U 10+ U 20+ U ~t+ “ 

s 

(3.23) 


After defining the various adjoint differential equations and boundary 
conditions by: 


8H 


X = on (to, 10). (10,20). (20, t ) 


8H 

* = -"x: on <W 

8x 


1;.U i =1,2,3 

l X 


(3.24) 


H (t f ) = 0 


equation for X „(t _) 
4 I 


X (t -) = X (t+) |f- ! + 4 

S S OX t - X 


X. (t+) ( i =1 ,5) 

l s 


s s 

H (t-) = H (t+) X (t-) 

8 8 OS 

Assuming expansion about a nonoptional initial control estimate, the 
quantities 5u and dm^ must be chosen to cause 6 J < 0. The particular 


6u dt 
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3Q_ lo A 

8*." 2 ? A J ' F" ■; 
1 V X 3 + X 4 


2 2 ' 
+X 5 


» i = 3 f 4. 5 


— = 0 

8x 6 


C. |T | =fcn(x 1 ) 


lili - 


A 

8x l 


Combining the developments above: 
Vertical Rise: 


dFj 0 T 8Q 9 C d 

8x^ ^ 8x, 


£s 

8x, 


= 0 


dF 1 8 Q 8C d 


8 *. 

t 


= 0 i = 4, 5, 6 


ifz 

8 x. 
i 


= 0 i = i 6 


8F, 

8x. 


= 0 i = 1 6 


(3.28) 


(3.29) 
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Pitch Over: 


9F x »T 


Cos C 


«« > -fe- vvl 

" K D M 


aF, 


9: 

9F, 

8x 


= 0 


»F, 

ax 


«f, 

!!i 

a*. 


i=-fse- c tQ !fel x 

4 K D 8 M 1 

i = .Its. c + q 5b1, 

5 K D ^sJ 5 


= 0 


8F. 


ST 


S ‘ **1 


Sin C, (t-t ) Sin C. 


vr 


Hi 

a*. 


= o 


8F r- 3 c 

2 = j8fi- c + Q 

5^ 1^3 D 8 *3.'"< 


3- 


3F_ 



8* 


= - Q C 


r? 3 - c. 


4 


D L 8X 4 


+ Q 


3C 

ax 


f] x < 


Hi 

3x c 


= Jaa. c to !lo] x 

|_8X S C D fQ 8 x 5 J *: 


C D*°>, 


!f?1 

8 «,J 


Hi 

3*, 


<3. 3C) 


(3. 31) 


0 
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0 








- o 
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Linear Tangent 

9F. 

— ~= 0 i =1,2 j =1 4 (3.36) 

j 


Second, consider the portions of the adjoint equations due to the 
performance index without the equations of motion adjoined (i.e. , due 
to the Lagrangian terms, L). 


A, Vertical - Rise, Pitch - Over, and Gravity Turn 

L =[q - AA(25)J 2 AA (37) 4 [aCC - AA(26)] 2 AA (38) 

where 

AA(25) = q max 

“■■"•I:* 

AA(26) = 3.05 

AA(38) =J P * [ACC-AA(260>O 
*0 [ACC -AA(26 jJ^O 
2 


4 h 2 -4 2 -s 2 ) 


q 

ACC = 


["V P<X l ,A .- qAC D] /x 6 

|^ = 2^-AA(25^ AA(37) + z[aCC - AA<26)| |£- CC 

8a. = i &-L 2 tx 2) 

fcCj 2 ax. r 3 + *4 + *5 / 


a* 

8 X - 


= 0 


(3.37) 


*2- - f r: 

3 


Ft a rr p - 

F-, •••'- i : 


" = fi x , 

•< *: , 4 


- AC T 


AA( 38) fc r (i-l. 



3 ACC 


= • [*» 3 Tj + A /x 6 

£C = r 8SL.1/ 

•4 L q »*, AC D 8 x 4 J ,x ( 


PACC 
dx 

8 ACC 
8 x _ 


.[a,^* , AC „ |a-]/ X6 

3 ACC 1 r . _ . . _ 1 

a = - — * I u I - PA - q AC _ I = 

8x 6 x | L sp e D J 


ftT 

5 ^ = - 2 [ ACC-AA ( 26 ) ] ACC AA(38) / x 6 
6 

B. Linear Tangent 

Lr = [ ACC-AA(27)1 2 AA(39> 

ACC = [ IT I ] /xt 
1 sp' 4 


8 ACC , ,**» 2 ,<v 

T=r — = - (u I ) /x = - ACC/x. 
ox, sp 4 4 


0 =1,2,3 

ox. 


■fk. = 2 f ACC-AA(27>] — ^ AA(39) 


= -2 [ ACC-AA(27)] AA(3?) ACC / \ 


Adjoint Equations for Vertical n vi rr * , Pi*~- 1 . O ver, 


ACC/x 6 


md Gra 
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Vfr + V -V'W 2 ! 

1 2 BF 

+ X 3 [ -vx 4 2 + x 5 V (Xj + R o ) 2 + 2k/ (x 1 + R o ) 3 + « 2 Sin x 2 4 ^/x^ ] 

+ \ [ - x 5 2 / f + r q ) 2 tan *2 1 f X 3 X 4 ^ ^*1 * R 0 ^ + ^ Sin x 2 Cos x 


^ /X 6 J 


(3.41) 


2 r ^ "I 9 p 1 

+ \ 5 [ x 3 x 5 /(X 1 + R o ) + x 4 x 5 / |^(x 1 + R o ) tan x^j + _J / x & ] 

-% z = ^- + X 3 [ 2 (x 1 + R q ) C? Sin x 2 Cos x 2 4 2 « x g Cos x 2 ) 

+ X - 4 [ “( x $ C so 2 x 2 /( x 1 +R o )+(x 1 +R o )«fcos 2 X;j - Sin 2 xJ-2»c 5 Sin xj 


+ X 5 [ * 4 x 5 Csc 2 * 2 / (x^ + R^) + 2 JJx 4 Sin - 2x 3 SI Cos x 2 ] (3.42) 

r 3F 

"*•3 = ^r 3 f S + x 3 [e^“ 1 x e] 

dF z 

+ V - X 4 /(x l + R o )+ i^ /x 6 ] 

3F 

+ X 5 [ -x 5 / (Xj + R o ) -2 S2Sinx 2 +9^ / x 6 1 (3.43) 

• a Y 1 

-V*r tx 2 'VV 

4 a F 

+ X 3 [(2 x 4 /(x 1 + R 0 ) +57 1 / * 6 ] 

4 

SF 

tX [ ' V!*! f V + SJ~ lx b 1 

4 

8F 3 

tx 5 [ - x 5 / [ (Xj T R 0 ) tan *.,] - ZO C / x 6 ] (J 44) 
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x 5 = i>r x 3 Ux s /<k i 4 v * ?QSin * 2 + ?T 7 / * 6 1 

5> *"> 

(• 

+ X 4 [2 x 5 /[(x 1 +R o )tan x^QCosX^-^p- /x fe ] 


8F. 


+ V - V i5 l + R o’ - V t <V V *“ * 2 1 ♦ sfA«(3- ' != -' 


-i. -X F/s 2 

V 6 ‘»x 6 S P l'6 

- X 4 




(3.46) 


Adjoint Egnatians for Linear Tangent Phase 

-*1 = + ^2 l ~~3 f ( *1 + R o }Z 4 2k,( *i * R 0 > J ] 

+t 3 : y 3 /(^ + v 2 i 

~ 8 ? + X 1 4 X 3 f ‘ *3^*1 f R o* " 

2 

+ ^> 2 [ 2 X 3 /^ + R q ) ] 

+T 3 i -J^/(^ + r o )1 

“V _ ik. C* S’ 2 'x ^ fZr ? - 

‘ X 4 ax. " X 2 F l' 5 ‘4 ' X 3 * 2 ' 4 
4 

3.7 Gradients 


(3.47) 

(3.48) 


(3.49) 

(3.50) 


The gradients with respect to- the control function, u(t), and control 
parameters c^, . . . , c,. are developed below. 


* 2 [ ACC-AA(26)1 

I AA(38)/W 

»F 

sp 

+ X 3 5T /x 6 

(Vertical rise, 

8F 2 , 

Pitch -Over, 

+x 4-sr '** 

Gravity Turn) 



r- ~r< 


II - 2 [ Ar C - AA( 27 ) ] 

u 

V X / ;; 

2 On ' A 

82T 

. ■ r 2 , 

+ X 3~ /x 4 


f i rTt>c?.r .* 




-X 


4 




The components for these gradients are as follows: 


(i) Vertical Rise: 


OF, 


1 * I 


On 


rn 


s h 

3" 

8F - 

¥T 


- o 


= o 




(ii) Pitch-Over: 


8F 

v — 1 = I Cos c ft-f I 

On sp 2 


ST =I sp Sin ^ tt - t vr ,SinC 5 

8F 3 

iT = 'sp Sin c 2 Cos c ? n. 


(iii) Gravity Turn: 


I'j 

On 


T. x„ 


/■ , t , 

*3 +x 4 ■‘•x. 


a Zi 


Byr 
" 3 


1 x i 

fip A 


4 


x„ + x . ■?*>:_ 

^ A u 


A X- 

-SR a 



50 


(iv) Linear Tangent; 


8F. I 
1 = sp 

9 u 




»F, I 
2 _ sp 

** TFT? 


c l : H c 1 = “ ^ x 60 + *6 (t o> " X 6 <*J> + X 4 ( V “ X 4 ^ 

8F »F »F 

c : H = X — — / x, + X r / x, + X c ~ — / x, 

2 3 8c^ 6 4 8 c 2 6 5 8c^ 6 


8 F. 


8c. 


~= -|T| [ Sin c 2 (t-t^fl (t-t vr ) 


8F. 


|T| t Cos c 2 (t-t yr )Sinc 5 


8 c. 


8 F, 

"icl 


|T| [ Cos c 2 (t-t„.)] (t-t_) Cos c 


vr' 


3‘ H = 


8 ^ 


l T i = aI sp - P< V A « 

- 8F 1 - 8F Z 

(X 2 /% 4> 8Cj + ('■ 3 / x 4 > 8 Cj 


sr* m n — v t 

K i (r 4 + 1 ) ,c 


8c, 


'4: H * 


!T| [(r ?I + iP /2 1 

8?, ^ 8F_ 

1 , ~ /X 2 

L / x . + X , r 

2 8c 4 3 8c, 

4 4 


/ ^ 
/x 4 


Sc , 


= | x | / (r 2 + 1) 3/2 


9“ * - |x| r/ cr 2 +i) 3/2 


(3.56) 
(3. 57) 


(3.58) 


(3.59) 


(3.60) 
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c 5 : 


H _ x jl !fi + x jl *Ii 

c 5~ x 6 *=5 + *6 Sc 5 


9F_ 

9c, 


& h 

9c, 


= 0 


9F 

8^ = |T| Sin c 2 (t-t vr ) Cos c 5 
5 


d h 

9 C C 


=-|T| Sin c (t-t ) Sin c t 


vr 


(3.61) 


3.8 Adjoint Function Boundary Conditions 


V 


x u =2 h (t f ^ ■ x u] p i 

X 2f = Z [*Z (t f ) ' X 2f] P 2 
Sf = 2 p3 {t f } ' X 3f] P 3 


H (t £ ) = 0 


4f 

‘s'' [ X 1 ^ 5 ] = [ X 1^2 X 3j 

r s T s 


9 gi 

wne re B. . = — . 

U 8xj 


1 0 
0 0 


0 0 0 
1 0 0 


(3.62) 


M x T (t s > 


B 31 B 32 B 33 B 34 B 35l 


H(tg) =H(t+) 


X 6 (t s ) 


(3.63) 


In order to calculate X. (t ) B ( j = 1, . . . , 5) and $ (t ) are required. 

r, 3 . x s 

j 

First, consider the equations for B^ ( j = 1, . . . 5). 
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Define: 


*1 1 *1 * R o 

g rJ x / +(* 5 t *i os ^ V* ={V 

-» — 


)M i 


'*-% -l ] 
-< ? ! 


2 2 ( 


) O Sin x_ 


B * = 


it a ( 


) Xj Q Cos 


B„ = ° 


B J4=I 


B 35 = l 


?[ f 

\Y z 
L .. L 


2 2 ( 


Finally, consider the contribution due to 6 (t-), where 
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COM Frr?:**.R *MZ>3URMKIU A X ~ON 
The computer tmplerr*er'‘-a fc *o-i of r> as»-NfK ton ?^ori f hms for 
somHo" cX Bolsa r-r'K'^m?: — -?.c •■ii^rrcsrf} h ritsr.re r ?. hi n»?pter 
3 the stntfte ascent opHsoir^HoD pmblem ms formulated and was not 

precisely in die Bolza form. la particular ;he performance index is 
composed of four integrals all of which ha m different equations of 
motion. Further, foe staging and filial times are not known bet implied 
by state variable terminal constraints, $ ( x^) = 0. Because terminal 
boundary conditions and state variable inequality constraints are bandied 
: - by penalty methods it has been found that a. considerable sa-ings in 
computer time can be achieved by real time human interaction with 
the executing program by way of a CRT display terminal. Problems 
' ; associated wifothe storage requirements of foe quasi -Newton algorithm--- 

also bad to be solved. One important clement of all the algorithms 

described in Chapter 2 is the one dimensional search, f !-D search). 

A large percentage of CPU time is consumed in the 1-D se? rch. Its 
»■ curacy and efficiency oave a great effect on foe success of the 
algorithms. These various practical considerations will be discussed 
in this chanter along wifo a full descriptors of foe nrog-a-n itself. 



trajectory opHr’i«ation prag^m. Sac smwbcc oi cherts is ssmnkr 

to those presented in Chapter 2 vr+h *+& addition of a CRT dispiay 
terminal for ho*o?n operator interaetif a with <he program vrhue it is 
executing. The main iteration loop consists of the forvrard integration, 
c backward integration, calculation of search direction. 1-D search, and 
convergence check. These operations are done repeatedly for a given 
set of penalty coefficients until an “acceptable" degree of convergence 
is obtained. At this point the human operator evaluates the final control 
and associated trajectory graphicaly on the CRT display. The penalty 

W 

_ \ 

coefficient values are appropriately changed by the operator and the 
iteration loop is re-entered. This process is repeated until the 
operator is satisfied that a further increase of penalty coefficients 
will sot yield a "better" control and execution is terminated. The final 
control and associated trajectory are then plotted by Calcomp for future 
analysis. 




OUTPUT 


Figure 4. 1 Flow Diagram 
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4. 2 Staging and Final Time 

The shuttle ascent trajectory optimisation problem is formulated 
in such a aty that t g and the staging and final times are free. The 
cutoff condition in both stages is determined by propellent exhaustion. 
Recall that 


Thrust = I u 
•P 

share u(t), the magnitude of the mass flour rate of propellent, is one 
of the controls to be optimised. All direct numerical methods require 
a first guess for u(t). This guests is usually stored pointwise at certain 
known times. Consider the boost phase of the trajectory. Assume u(t) 
is stored at n equally spaced storage locations and is piecewise linear 
between storage locations. 



Figure 4. 2 Control Storage 

Uj (i = 1, n) and n are known values. 

The mass of propellent is known and must equal the area 
under the u(t) curve, 

v[r'V' • • + Vi + -jp-] “ 
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Thus 


m„ 


[H "t.1 

[— + “2 + Vi T J 

Since At is known, t can be computed, 

8 

t = (n-l) At (assume t =0) 
s o 

The same procedure is used to calculate t^. 

Now consider the interval between storage locations J and J + 1 

which corresponds to the time interval tj to t j + 
u 





Figure 4. 3 Single Control Segment 

On this interval, u(t) = 2 a T + b where T = t-tj and 


b = u. 


2a - 


u » a 

j + i j 
At 


Since u = -m, the mass can be obtained as a function of T by integration, 

T 

m (T) = m_ - f ( 2 as + b) ds 

J k A 


* A T +BT + C 


Where, 

A = -a 
B = -b 

C = m ( mass at start of interval) 

J 

Because of the assumed form for u (t) it is possible to calculate 
t^, t^ At, and mass (t) analytically. This avoids the problem of guessing 
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a tg and At and integrating the whole system of equations while checking 
for fuel exhaustion. It also avoids the need for extending or contracting 
the control guess if the mass of propellent is not zero at the guessed tg. 

4. 3 Storage Problems With Cuasi-Newton Algorithms 

h Section 2. 2 it was shown that 2 i + 4 time functions mus t be 
stored after the i** 1 iterate in order to compute the i + 1 st search 
direction. Each of these functions is stored as a n- vector of numbers 
which correspond to the function values at n equally spaced points on 

[t fl , tg ] . Thus ( 2 i + 4) n floating point numbers must be stored after 
th 

the i iterate. The computation per iterate also increases because 
of the increased number of inner product evaluations. Thus it is a practical 
necessity to restart the algorithms to a pare gradient step every iterate, 
ft has been found that 3 < q < 8 is a good choice. The value of n must 
be larg. enough so that a'feood" representation of the functions is 
obtained. For the shuttle optimization problem the time interval is 
approximately 500 seconds and n was chosen to be 500. Thus storage 
must be allocated for (2q+4)n=(2x8 + 4) 500 = 10, 000 double 
precision floating point numbers. Additional storage must be allocated 
for other variables used in the program and for the object program which 
is generated from a fortran source aeck of 3800 statements. 

During the initial testing of the program on the University of 
Michigan IBM 360/67 virtual memory computer all storage was done 
in fast memory. Hovrover it was found that core storage was exceeded 
when the program was first run on the JSC' s Univac 1108 ■ .nputer. 



60 


To overcome this difficulty file 10, 000 double precision floating point 
numbers needed for the quasi-Newton algorithms were placed on drum 
storage. This reduced the amount of core storage required allowing 
file program to fit on the 1108. Upon running the modified program on 
file IBM computer a considerable savings was realised in reduced 
v irtual memory charges, ft was also found that no significant increase 
in file amount of CPU time was incurred. There are two reasons for this, 

i) a very small percent of CPU time is spent calculating file search 
direction. Most of the CPU time is spent integrating the equations of 
motion. On each iterate a forward integration and a backward integration 
are required to determine the gradient and a number of cost evaluations 
also requiring forward integrations are performed by the 1-D search. 

ii) the updating equation for H. y. and the equation for d. are 
summation . j which require inner products of the stored functions in 

the Race s^qumce as they were generated and stored. Assume H. ^ y. ^ 

and d. s.re to he calculated. H y through H. _ y . _ are stored in a file, 
i o o x -2 i-Z 

Start of 
File 


/ 


0 

a 0 

EM 

v 2 

1 

1 

1 

H i-2 y i-2 


T 


\ 

Read 

Write 

Pointer 


Figure 4. 4 File Storage Diagram. 
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At the end of the last iteration the file has been rewound. The updating 

equations for K. , y. , will read H y , H,y., . . . , H. * y . , in order, 

1-1 i-l o o 11 i-Z i-Z 

calculate H. , y. ,, then write H. . y. . onto the file and rewind. Con- 

currently the equation for d. has been using the Hy functions. The files 

in which Hy and s are stored need only be rewound once on a given 

iteration and no forward or back spacing is required. Even if tape 

were to be used as the storage medium, instead of fast core storage, 

the increase in computer time would be small. When drum storage is 

used the increase in computer time is insignificant. Thus there is no 

need to restart to a gradient step because of limited storage. 

As mentioned previously the computation time per iterate increases due 

to the increasing number of inner product evaluations which must be 

made. The inner product is a quadrature. 

/>*£ T 

<u,v>= J u v dt 
t 

o 

where a and v are stored pointwise. If it is assumed that the stored functions 
are linear between storage locations die evaluation of the inner product 
reduces to a summation. Consider the interval t^ to t^, 

v 

*1 *2 

Let T = t-tj and Lt = t_ -t^ then on J 

u (T) = a T + b and v (T) = « T + p 
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where 


a = V!i 

At 


* ~ V 2 " V 1 


At 


b = 


P = v i 


The inner product of the functions between ^ and ^ is; 


<u, v: 


f At 

> = / (aT + b) (<r T + p ) 

\ z 0 


dt 


At 

= / faff T 2 + ( <* p +*b) T + bp ] dt 
0 


= **r At 


At + bp At 


and die total inner product is 

n-1 

».▼>*,.* 53 <u,v> 

‘_l W - _A t.*I . , 1 

o f 1=0 i l+l 


<u. 


It was found that this method of evaluating inner products is considerably 

faster then higher order quadrature formulas and that convergence rates of 

tta algorithms do not suffer. 

4,4 One Dimensional Search ( 1-D Search) 

On each iteration a search direction d. is generated, and then a 

new control is calculated, 

u . . . = u. +Cl. d. 

i+1 l it 

The goal of the 1-D search is to find a scalar parameter a ^ which yields 
the greatest cost decrease. At such a value it is necessary that 

J ( u. + ff d.) = 0, 

da l i' 

which is an important element in the convergence proofs of the quasi - 
Newton algorithms for the linear quad, at ic problem. 



63 


A large fraction of CPU time is spent within the 1-D search 

and its accuracy and efficiency greatly influence the convergence 

a J 

rate. For small cr, a a < 0 thus we can expect a functional relationship 

with the following form. 



Figure 1. 5 Cost vs or 

As or increases J will decrease until the higher order terms in 

the expansion dominate and J begins to increase again. The 1-Dsearch 

attempts to find or *. Th performance index is evaluated at or = 0 and 

* 

or where or. is a guess for a . a is then increased or decreased 

until the minimum is bracketed, that is three points are found sueh that. 


J(«> 



cr^ ^ a 


J(or . ) > J(ff •) < J ( a t) 
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The function J(« ) is approximated by a quadratic curve, 

2 

J (a ) = a or + b or + c 

where a,b, and c are determined by fitting the quadratic curve through 
the three data points. The minimum of the quadratic curve is given by, 

~ = daa + b = 0 
8a 

-b 

* = 2a 

The performance index is evaluated at or 1 and if, 

J(a') < J(* } 

■r 

the control generated by or is chosen as the local minimizing element . If 
J(«*) > J( ttj ) 

a new quadratic curve fit is performed with a' , a ^ and, 

a . if a ' > a . 

1 3 

or 

a, if a' < a . 
k J 

This process is repeated until a minimum is found. 

4. 5 CRT Graphic Display 

The. space shuttle ascent trajectory optimization problem developed 

in Chapter 3 is a Bolza problem with the addition of state variable 

inequality constraints, 

acceleration < ACC 

— max 



Dynamic Pressure < Q 

~ max 

and terminal boundary conditions, 
50 x 100 nm orbit 
Inclined 28.5° to equator 
Entered at perigee 


whe re 


ACC = 3. 05 g' s boost phas.. 
max & 

ACC = 3.0 g's orbiter phase 

max 

Q = 650 psf. 

max 

This optimization problem is replaced by an unconstrained optimization 
problem where the terminal boundary conditions and state variable 
inequality constraints a e enforced by the method of penalty functions. 
The . aw unconstrained optimization problem has seven independent 
penalty coefficients, and '■he performance index is 


J = - W^ + Pj [ ERROR IN FINAL RADIUS]. 2 

+ P [ ERROR IN FINAL RADIAL VELOCITY ] 2 
+ P 3 [ ERROR IN FINAL TANGENTIAL VELOCITY] 2 

+ P 4 / s (q(t) - Va/ u W'-W d ‘ 

o 

+ P 5 / s <acc(t)-acc ma / U ( *cc(t) - d» 

*o 

r t f 2 

+ P, | (ac c(t) -acc ) U (acc(t) -acc ) dt 
6 J max max 

s 

+ P ? [ ERROR IN FINAL INCLINATION] 



Here P. ( i = 1,2, 3, 7) are penalty coefficients associated with die terminal 
boundary conditions and P. ( i =4, 5,6) are penalty coefficients associated 
with the state variable inequality constraints. For a given set of penalty 
coefficients a particular uncor strained optimisation problem is defined. 

The solution to the original constrained optimisation problem is 
appro xim a te d by a sequence of solutions to the unconstrained problem 
generated by letting F. ( i = 1, . . . . , 7)-*w. As P ( i *1,2, 3, 7) are increased 
the solutions generated will more closely satisfy the requirements of a 
50 x 100 nm orbit inclined 28. 5° to the equator en.ered at perigee. 

Likewise as P. ( i = 4, 5. 6) are increased the state variable inequality 
constraints on dynamic pressure and acceleration are more strictly 
enforced- The ultimate goal is to find fire control history which yields 
the maximum liftoff weight and satisfies all seven of the constraints. As 
expected, in practice as one oenalty coefficient is increased the error 
associated with it will decrease while the errors associated with the 
other coefficients will increase. Thus by improving the trajectory 
in one respect it is possible to lose something somewhere else. 

-Sensitivity to changes in fire different i?«t* *ty coefficients also varies. As 
the penalty coefficients become larger the overall problem will become 
increasingly sensitive to changes in the control and numerical instability 
will eventually, resirlt. The way in which the penalty coefficients are 



increased will strongly influence the overall convergence rate of tfte 
algorithms. The main drawback to the method of penalty functions is 
that the penalty coefficients must be increased in a problem dependent 
way. Even for simple example problems which require little computer 
time for a trajectory integration and which have only one or two penalty 
coefficients, the choice of these coefficients and the way in which they are 
increas«*d is critical for rapid convergence. Because of the complexity 
and relatively long compute r time required for a trajectory integration 
of the shuttle ascent optimization problem a better method than trial 
and error is required for chooahg penalty coefficient values. 

By using time sharing computers and CRT display terminals the 
problem of choosing penalty coefficient values can be very efficiently 
solved by human operator interaction with the executing program. At 
the end of eacS iteration execution is terminated and control transfered 
to a CRT d i splay terminal. Because of time sharing this interruption 
of the executing program is very inexpensive. At the request of the 
human operator important information is then graphically displayed on 
the CRT. The information is evaluated and a decision on changes of the 
penalty coefficients is reached. This information is communicated to 
the computer and execution proceeds. By placing a human operator in 
the program iteration cycle convergence times are red ced, the 
computer is used more efficiently, and the operator quickly builds an 
intuitive feel for the physical problem being solved. 



For tile shuttle ascent optimisation problem it is helpful to 
graphically display dynamic pressure, acceleration, and u as functions 
of tune along with terminal miss values. The best convergence rate 
eras achieved by first increasing P. ( i = 1,2, 3,7) yielding a trajectory 
which comes" close' 1 to the desired terminal boundary conditions . Then 
P. ( i • 4, 5, 6) are increased to enforce the state variable inequality 
constraints while simultaneously increasing P. (i* 1,2, 3,7) so that all 
intermediate trajectories remain "dose" to the terminal boundary 
conditions. 

The anility to interact with the executing program can be useful 
in other ways. The interrelationship of adjoint, state variable, search 
direction., and gradient time histories can be conveniently a n a l y s ed using 
the CRT display, in conclusion the ability to communicate with the 
executing program is a valuable tool for analysis of optimization programs. 

4.6. Subroutine Description 

The computer pregram consists of nineteen subroutines controlled by 
the main control program. Figure 4.6 presents a subroutine map which 
illustrates the relationship between subroutines, hi this section the 
function of each subroutine will be explained. 

MAIN - reads input parameters, calls SPLINE to obtain curve fit of 
aerodynamic coefficients, controls forward, backward, and 
cost integrations, cadis CAL to determine constant gradients, 
calls SEARCH which contains the 1-D Search, also contains 
logic for interaction with CRT display terminal. 
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* tho re are three FCT' 8 and three OUTP' 8 I = 1, 2, 3 


Figure 4. 6 Subroutine Map 
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SPLINE 

INTEG - 

ALGOR - 

CAL 

SEARCH 

TRUNC - 

LAMF 

LAMS 

POLAR 

DRKGS 

FCT 

OUTP 

CRAFT 


• c ub routine which interpolates by piecewise cubic splines 
aerodynamic coefficients such as which is given in 
tabular form as a function of Mach number. 

contains the logic to determine the mass distribution, staging 
time, final time, and caiis DRKGS for forward, backward, 
and cost integrations. 

contains the various algorithms which require a gradient 
g{t) as the input and producb d(t) the search direction as the 
output. 

performs the quadrature which calculates the constant gradients . 

- contains the 1-D search, i.e. determines a which minimises 
the per for mance index, see Section 4.4. 

performs the truncation of new controls generated by varying 
a in SEARCH. 

• calculates foe value of the adjoint variables at t^. 

- calculates the jump in adjoint variables at t^. 

• calculates foe jump in state variable at t and calculates foe 
inclination of foe orbit. 

. a double precision fourth order variable step sise Range -Kutta 
integration subroutine contained in foe IBM SSP package. 

- computes the right hand side of the system of equations to 
be integrated. 

- an output subroutine used by DRKGS 

dC 

- calls spline to determine and P • 
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ATMOS - calls MODEL to determine density (a), pres sure (p), and . 
speed of sound (a); also calculates 
dP dp da 

dh ah dh 

where h is the altitude. 

MODEL - contains the atmospheric model; see Appendix B. 

4. 7 Numerical Results 

The final control history and associated trajectory which will 
be presented in fins section are the result of computer runs made at 
both JSC and at the University of Michigan. The initial control guess 

Cj = payload mass = 80, 000 lbm. 

C 2 =7 = 0. 689241° /sec. 

C 3 = a = -0. 431410 x 10" 3 

C 4 - b = 0. 365070 

Cg = i|i = -19.0049° 

u (t) = 98? 

Tbi s resulted in a trajectory with the following terminal miss 
values, 

A R = -180, 000 ft. 

A U * -200 fps. 

A V = 602.1 fps. 

Inclination = 26.23° 
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The staging time was 118.73 sec. with a final time of 503. 3 sec. 
The state variable inequality constraints were also violated. Q reached 
a peak value of 819 psf at 65.2 sec. while the maximum acceleration 
during first stage was 3. 8 g 1 s and during second stage was 3.9 g' s. 

With the above initial control the program was run on JSC' s 
computer for 12.75 minutes. The resulting final control became the 
initial control for subsequent runs made on the University of Michigan 
compute'’ An additional 45. 7 minutes of computer time were expended 
on the University of Michigan computer for a total run time of 58. 4 
r v.oses. The penalty coefficients were; 



INITIAL. 

FINAL 

P 1 

10 6 

IO 16 

P 2 

io 9 

10 W 

P J 

io’ 

10 W 

P 4 

10* 

io“ 

V 

10° 

10 3 

P 6 

IO 10 

1<P 

P 7 

io» 

io 15 


The final control is; 

C = 101, 300 lbm 
C 2 = 0.631857 °/sec. 



C, = -.47854x10 
C 4 = 0. 366590 
C 5 = -8.5° 
u (t) - Figure 4. 7 

On the converged trajectory, the staging time is 121.1 sec. and the 
final time is 504. 0 sec. Figure 4. 8 shows the angle y above toe local 
horisontal at which the thrust is orientated. Figures 4. 9 and 4. 10 indicate 
that the state variable inequality constraints are being enforced. Figure 
4.11 shows the time history of altitude vs time. The terminal miss 
values are, 

ft R = -4,700 ft. 
ft U -1.2 fps 
ft V- 5.2 tos 

<r » 

Inclination = 28. 8° 

These values could be improved by decreasing toe integration 
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CHAPTER 5 


TERMINAL CONS "RAINTS AND VARIABLE t INAL TIME CONSIDERATIONS 

When a function space gradient-type technique is employe ), me 

must continually confront the problems associated with terminal ■ otistr wV- 

and variable final time. Usually for flexibilit* .nd ease of progranw. ; og 

a penalty function type of approach is preferred to a projected gradient 

approach. In this chapter we shall compare the standard penalty function 

- 1ft 

ap r-ach to the Hestenec multiplier method on a single problem. The 
results for this particular problem indicate that the multiplier method 
does not improve the performance enough to merit the additional programmer 
and complexity. 

The major result of this chapter is concerned with a simple procedure 
which apparent!^ improves considerably" the rate of convergence 
gradient-type methods v hen penalty functions are employed on -.nibble 
final time problems. The result is simply that the initial estim .re of 

(o) . ' - . - * : 

t » say t ' , should be less than the optimal t , value,, say t . 

f f r _ a 

fa fact, it appears that a physically unrea conable choice of 

(01 * - ’ 

which guarantees t^ < l f is superior to a physically reasonable ini'ial 
( 0 ) ^ - 

trajectory with t f > with respect to rate of con ver gene Although 

this property has yet to be proved mathematically, it apy- - • '.-.be 
heuristically justifiable, and all rf our -umerjcal simulation, confirm 
the trend. Finally, it will be shown that a recently pror >sed method for 
treating variable final time problems by Tripathi and Narendr,.^ is 

"0 
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essentially th~ well -known method developed by Long *' in 1965. 

5.1 Variable Final Time Prob 1 ;ms 

Consider the performance index for a time -optima’ control problem 

with fixed -endp^ «t type terminal conditions 
i* n 

- T ■ c *i *5i p i "i - ‘u' 2 + E c i <x i - v (5 - 1 » 

where 


P. are penalty coeffic len.s for terminal constraints; P. = 0 if x. (t,) 
i t i ' r 

is not specified. 

18 

C. are multiplier constants for the multiplier method ; C. = 0 if 
♦he penalty function method is used. 

In the time -optimal control problem die algorithms require an 

CO) 

initial estimate of t^, say t^. On future iterate?, a procedure for 
updating t^ must be specified, and this will affect the rate of convergence. 

The following have been proposed in the literature. 

1. If t. is reduced, the pertinent functions are well-defined for all t. 

If increases, the control i? set equal to the value at which 

is of the previous iterate, in the extended interval f tjj^ t, . 

The program used in this chapter is based on this technique 
and it worked satisfactorily, at least for relatively 

simple pro v '"rms. 

2. If increases, the various functions are suitably extrapolated 
ov~r the new range. 

3. The functions maintain same form, only tlie time scale is 
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modified to tike care of the changes in the interval f t^_, t^]. 
Tnpathi and Narendra ^ found this method to be satisfactory 
in practice. 

4. Convert the problem to a fixed final time problem with the 
transformation below. ( An additional parameter appears due 

17 

to this transformation, and the method was proposed by R. Long ). 

A. Define t = a s 

where x = f (r, u, t): the equations of motion 
a = constant to be determined 
s = a new independ int variable. 0 _< s _< i 

B. Let: s = 0 be the initial point. 

s = 1 be the final point = a. 

C. New equations of motion 


( ) 



x = af (5.2) 

a' =0 

Note: Since "a" is an unknown constant parameter, 
an initial estimate of ,: a" is needed and an 
updating scheme must be defined. 

It will be shown that die r*ethod proposed by Tripathi and Narendra 
is essentially the same as the method by Long. The justification is as 
follows: 


If t 


(i +D 
f 



then following Ref. 16, a function. 


say k {•:), 


ah' aid be updated by 
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C l+1 >(t) = k (l) (t/ft). 


( 5 . 3 ) 


This method causes the function k (t) to be compressed or 
expanded as t^**^ d< creases or inc~sases, respectively. 

The method can be defined alternatively by introducing the independent 


variable T - a t, Then, the function for the next iteration is 

k (l+1) (t) = k (l) {t/a ). (t =0) 

o 


( 5 . 4 ) 


For example, if 

k l (t) =t 2 +t 


V t« 


[° • '*] ’ 


then, the function for the ( i + 1) - iterate will be 

2 

k ‘ + 1 (t) = k 1 (t/a ) =^- V t * 0 , a t^ . 

Thus, Tripathi <md Narendra' s method can be represented by the 
transformation T = a t, if the relation t^ * + *= a t^ 1 holds. 

In the application of Long* 8 method, the value of the constant "a" 
has to be guessed initially to start the scheme. Assume that at the (i +■ 1) 
iteration. 


i +1 i 
a = a a 


(5-5) 


This implies 


i + 1 i + 1 
« = a s 


(5.6) 


= a s. 


(5.7) 


Substitution of (5 -*to (5.6), and use of (£.7) implies 

i + 1 i + 1 i it 1 i 

t = a s = u a s = a a — r = a t , 


„ i +1 i 

t = a t 


( 5 . 8 ) 



which is the transformation equation used by Tripathi and Narendra. 

Thus, method (3) and method (4) have similar basic characteristics. 

hi the n.jtt section meAod (1) above is employed, and numerical 

( 0 ) * 

examples are presented to show Cat t £ < t £ gives a more rapid rate 

( 0 ) * 

of convergence than t £ >t £ . 

5. 2 Numerical Examples for Minimum Final Time Problems 
Example 1. Zermelo 's problem 
x = v cos 0 
y = v sin 0 

1 

e =u (5. 9’# 

where v = constant, 

x(0)=x o = 0. y(0)=y o = 0. 0 (o) = B q = 0 
|u | < k, k maximum turning rate 
Determine the minimum dme to reach the specified final states: 


A. 

x (t f ) = free, y (t f ) = free, 0 (t f ) = 0 £ 

(5.10) 

E. 

x (t £ ) = x r y (t £ ) » y £ , 0 (t £ ) = free 

(5-11) 


For these simple problems, analytical solutions can be obtained 
without difficulty. 

Case A: It is easily shown that given v = 1, k = . 50 and 0 £ = 2v, the 
control will be eithex u = + k or u = -k for the vehicle ‘r rea- v the 
specified heading in minimum time. 


The cost functional is 
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J = Ct f tP^X-Xj) +P 2 (y-y £ ) +P 3 (e-e f ) 


+ C £ (x-x £ ) + C 2 (y- y£ ) 4 C 3 (0-0^ 


(5.12) 


where P £ = P 2 = 0 and = C 3 = 0 for the penalty function method. 

The optimal trajectory is a circle centered at (0, 2) with radius two for 
a a 

0 £ = 2w = . 5 t £ , t £ = 12. 56 second?, and 


Initial Conditions 


Terminal Conditions 


x = 0 
o 

*. = 0 

0=0 

o 


x £ = free 
y £ = free 
0 £ = 2x = 6.28 


(ii) Consider t £ ^ = 19 seconds ~>> t £ . After six iterations, l^= 


The purpose of this problem is to show how the conjugate gradient 

( 0 ) 

method is affected by the initial final time estimate, t £ . Let C = 1, 

C 3 =0, Pj = 100, and the integration stepsixe = At = 0.2 seconds. 

(i) Consider =2 seconds « t, . The algorithm increases 
f a f 

(21 

the final time to 12.65 seconds with u = -K 5 in two iterations, 

(see Figure 5.1.a). 

a /A i 

12.084 
/2\ 

seconds with u = +.5. (See Fig. 5.1.b. ) After two iterations, t^' '= 18.96. 

(0) * 

Thus, both cases converge lapidly, with the ^ = 2 ' t £ case 

having Ae fastest rate of convergence. 

Case Be The exact solution for this case is as follows : To reach the 
specified position* in minimum time, the vehicle will first turn at the 
maximum rate, and then switch to the singular arcu = 0 for straight 
line flight to the desired position, i.e. , 

u = + k Y t « f fc 0 . *i J 

u = 0 Y i * {tj, t“| (5.13) 
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u 


.5 

.4 

. 3 
.2 
1 

0 


\ 

\ 

\ 

\ 

\ 

\ 

\ 

\ 

I 

I 


4 £ 


o 

— initial control estimate, tf =2 

optimal solution 

control at 2nd iteration 


& 10 12 14 U 20 22 


*■ 


(a) 


f (0) 

f 


a 

f 



control at eig^fc iteration 


2 4 6 8 10 12 14 16 18 2022 


(b) t f (0) » t* 


Figure 5.1 Control Profiles for Example 1. Case A. 


The • ■orrespondipg cost functional is (5.12) but with C = 1. 


P. = P - Mf;, - C - P = C = 0 (penalty function method) 

i c l t* 3 j 

Initial '.onditions Terminal Conditions 


x = 0 

n 

v =0 


= 0 


= 4 miles 
y = 3 miles 
0^ - free 


The optimal final time is t* = 5.058 seconds. Four cases are 
considered in this example. 

(i) t^^= 2 seconds <x t * . The final times goes to t^ = 4.95 
seconds on the first iteration, and to tj^^ = 5. 042 seconds 
in ten iterations. The position error after ten iterations is 
within one percent (see Figs. 5. 2. a and 5,3). 

(ii) = 4 seconds. This guess is cion.* to the true minimum. 

The program performs smoothly and the terminal position 
error is less than one percent (see Fig. 5.2.c). 

(iii) t^^ = 6 seconds ( slightly '„rger tuaa the true minimum). 

( 12 ' 

Little improvement in final time, t =5.78 seconds, after 
twelve iterations. The position error is about 2.5 percent, 
and the program terminated due to insignificant cost change. 
Another interesting aspect of this case is that the control 
profile converges to a profile far from, the optimum. This 
implies that an initial guess with t^^ > t* may have the tendency 
converge (apparently) to nonoptimal solutions (see Fi^s. 5.2.d 
. ..3). 

(!) * 

► =10 seconds » . After eight iterations, the position 

iQ\ 

re - 1 is less than .2 percent. However t ' =9.87 and again 




initial guested trajectory t / 10 


! 

1 



the control profile moves in the wrong direction (see Figs. 5.2. b and 5. 3). 


Example 2 -Flight in a Horizontal Plane 

The equations of motion for a coordinated turn in a horizontal 
plane, with the thrus* always aligned with the velocity, are 
dx 

— - = V cos b 
dt K 


£ * v i» 

dV (T-D) 


dt m 

v dp _ L sin a 

dt m ( 3 . 14) 

dm _ -X 
dt c 

To maintain the aircraft in the horizontal plane, an algebraic constraint 
is imposed 

L cos <r = mg, (5. 15) 

and the parabalic drag polar is assumed, i.e. , 


C - C + K C 
D DO L ’ 


(5.16) 


where and K are independent of the Mach number and the Reynolds 

number and 

L = SC T V 2 


d4p SC V 2 . 
2 K D 


(5.17) 


Two of the three functions T,fl ,<r may be identified as controls with 
the Lhird determined by the couscr^nt (5 15). In this example thrust 
magnitude and bank angle are controls which are all bounded, i.e.. 


•)0 


(t) . < c r >'t) < r (t) 

1 y mir, — — max 


t . < r < i 

min — — max 


(5.18) 


a (t) . = -nr (t) for a symmetrical aircraft, 

mm max 

The cost functional to He minimized is 

J = Cr / 1 P 1 (X ' X / + P 2 (y ' y f } 2 + P 3 (V - V f )2 + P 4 < P -P/ + P 5 'm-m f ) 2 

+ Cj (x-x f ) + C 2 (y-y f ) + C 3 (V-V f ) f C 4 (p -p f ) - C ? 

(5.19) 

A relatively simple case is selected to show how the initial final 
time estimate will affect the performance. 

Initial Conditions Terminal Conditions 


x =0 
o 


' o 


- 0 


V =2.2 Mach -^2136. 2 ft/ sec 

a - 

o 


= 6 miles 

= free 

= free 
0 - free 


W =■ 434 lbs 
= 1000 , 


w = 861 lbs 

Again the penalty function method is used, and = lOOO f P^ 
C. = 0, i = 1, . . . , 5, C = 1. 


P 3 * P 4 * P 5 * °- 


The optimal solution for this case is t^ = 5.7 5 seconds. The thrust 
profile is the boost-coc.st f ype and the bank angle is zero for all time. 

Three final time estimates are considered. {Figure 5) 

1. t - 4 seconds, program forces the final time to the iborhood 

of t* in three iterations and obtains = 5.79 on the fourteenth 
iteration, with f he boundary condition error less than 0.2 nercen 1 . 


(see Fig. 5. 4. a). 



V 


( 0 ) 


6 seconds. This h auain in estimate close to optimal, 


hut slightly longer than the true minimum. The control profile 

(14) 

approached the boost-coast type and ended up with - 5.94 


in fourteen iterations, ^ee F’g. 5. 4 .b). 

/ rt i 

3. ; - 10 seconds. After twenty - three iterations the terminal 

position error was less than .1 percent (x^. = 6.00003), hut there 
was insignificant improvement in flight time, t^ = 9*85. At the 
thirtieth iteration, started to improve to 8.76 which is still 
far from the true minimum time. The thrust control profile 
tends to the coast-boost type which is far from the optimal 
solution (see Fig. 5. 4. c). Heuristic reasons for the behavior 
in the examples above are given in Section 5.3. 

5. 3 Method of Multipliers 

A brief comparison of the penalty function method and the modified 

multiplier method (M M - 2), Ref. 19, was undertaken in the study. Based 

18 

upon the theory by Hestenes , M M - 2 should ^ jrform better than the 
penalty function method. Our experience has been that, with the conjugate- 
gradient algorithm, some improvement does occur. However, the 
improvement is not significant enough tr justify the additional programming. 
5.4 Conclusions 

The examples in Section 5.2 demonstrate numerically that the 
(Q) 

choice t^ v appears to improve considerably the performance of 

gradient -type methods when penalty functions are employed. Although 



2 





rigure 5.4 Confrol profiles for Example 2. 




ix mathematical proof of this fact (which would involve- a rate -of - 
convergence-type proof) has not been obtained to date, the fo ./mg 
heuristic argun ent is offered in suuporc of- the' po$sltle7*ieneraUty o r 
this observation, 

Con r .ider a time optimal control problem V it h t' v il constraints 
where t* is the optimal final time. Suppose . -Theii, -t is im- 

possible for the initial trajectory to meet, tho bovmdai conditions', and 

oust be greater than to decrease the error on the terminal 
f i . 

constraints. Th'*s> the optimal solution has the unique characteristic 

of being the closest trajectory to the initial iterate, with respect to 

final time, which satisfies the terminal constraints. On the other hand, 

if > tjp , then it is probable that there exist infinitely many nearby 

solutions which satisfy the terminal constraints. Sin. ° with penalty 

functions terminal constraint satisfaction is a major part cf the per- 

fo nuance index, thera exists the tendency to r 'lock-ir on the terminal 

conditions at t > That is, the optimal sohition m, longer possesses 

the unique property of being tie closest trajectory v'V^h satisfies the 

boundary conditions. With regard to mathematical imolications , the 

(0) * 

statements above imply that the mini m is ’flatter 1 ' 'ft > t than 



CHAPTER 6 


THE PRAXIS ALGORITHM 


fin the previous chapters function space algorithms for minimisa- 
tion have been studied, fin this chapter we shall consider a recently 
developed parameter optimization scheme which does not require the 
objective function tc be differentiable. Such a scheme is of use in 
problems where it is difficult or even impossible to find the partial 
derivatives of tbo objective function directly. 

20 

\ e shall 'il.st discuss Powell' s method and the modifications 
23 25 

due to * letcher and Brent . Some specific- properties which are 
closely related to convergence are presented along with an application 
of the method to a time -optimal control problem. Also, the subroutine 
of Appendix D has been built into the NASA- JSC PEACE parameter 
optimization program. 

6.1 Powell* s Algorithm 

The basic concept of Powell' s Algorithm is to minimize a scalar 
function of n variables, say f (x* , .... x & ), by searching along n 
directions which span the space. Thus, for one iteration, the basic 


procedure is as follows: 

Let Xj be the estimate of the vector x on the iterate, and ty 
be vectors which span the space ( initially [ u„] is the nxn identity 
matrix). Then: 
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1. For i - 1, .... n, compute 0. to minimise f ( x_ . . f 0. u.l 4 

l Ji 1 ™X & \ 

where x_ . * x_ . , + 0. u. end x_ _ S x_ 

J,t J, l-l ii J, 0 J. 

2. For i = 1, . . . , n-1, replace u. by u. , 

i i r l 


3. Replace u by x„ - x_. 

a J, a J 

4. Compute 0 to minimise f ( x_ + 0e ), define x _ . . * x_ +<*u 

J n J r l J i 

return to (1). 


A simple graphical example will clarify the iteration procedure. 


Consider an ellipse in two-dimensional space (see Fig. 6.1). The 


algorithm starts at ( x , y ), and on the first iteration, the search 

o o 

directions , u^ are (1,0) and (0,1), which are along the x and y directions. 



Fig. 6.1 Operation of Powell' > Method. 
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Following Step 1, the algorithm searches along the 
direction, and B is the resultant point along the - axis. Then, 

the algorithm searches along the bl ^ direction from the point (x + d u/^\ y ), 

6 O 1 l O 

and C is the resultant point. Steps 2 and 3 require the new search 
directions to be: u^ 1 * = u 2 * 0) , = (/> \if°\ where u^* 

is actually the direction along AC. Finally, Step 4 implies that point 
D is the value of the new estimate, i. e. , x^ = x q + /iu^ . The procedure 
is then repeated, and for this example, convergence will be obtained 
an the second iterate. 


Ob the second iterate, the new direction O D* is conjugate to u 
according to a theorem developed by Powell^. The minimum is obtained 
by performing an additional search along this direction. 

This example demonstrates that the algorithm converges in a 


finite number of iterates for quadratic functions. As one might expect, 

the property of conjugacy plays an important role in this connection, 

and more details will be presented in the next section. 

This basic procedure has the defect, as pointed out by Zangwill 

that a poor guess of the initial position (e. g. , point B in Fig. 6.1) might 

lead the algorithm to fail to find the minimum. Instead, the algorithm 

will converge to a minimum along the line u^**' . which defines a proper 

2 

subspace of the space R . 

20 22 

hi order to overcome this difficulty, both Powell and Zangwill 
proposed methods to retain the linear independence. Numerical 
experiments in Ref. 25 show that Powell' s modification is preferable 



to Zangwill' s. 

6.2 The Roles of Conjugacy, Orthogonality and Independence 

By definition , Nro vectors and ^ are said to be conjugate 
with respect to the positive definite symmetric matrix A if 
«1 T A <*2 = 0 

A set of conjugate directions is a set in which the vectors are pairwise 


conjugate. 

2 2 

Consider a quadratic function f (x, y) = * + 4y , with elliptical 
candours as shown in Figure 6.2. Writing the function in matrix form, 


* (*.▼) s l *.y ] [: :i;l 

(4 


with the A matrj 
A 


and the two unit vectors along the x and y axes 
*1 


two unit vectors al 

■tl •■■■[:] 


(6.3) 


(6.4) 


(6.5) 


It is obvious that and u^ are conjugate. 

With Powell' s method one can obtain the minimum by searching 

along each of the conjugate directions once as long as the space on which 

the function is defined is spanned by the conjugate vectors. From Figure 

6.2, one can easily see that the minimum is obtained by searching along 

the x- and y - directions only once regardless of the initial guessed position. 

o 

Now consider the same function in a coordinate system rotated 45 


from the original system: 



Figure 6.3. Convergence characteristics for nonconjugate and 
and conjugate search directions. 
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* (*. y) = [ *. y ] 


1 - 2 

2*2 

2 2 

2 2 


( 6 . 6 ) 


From Figure 6. 3, successive searches along the x and y axes 
will not reach »he minimum due the fact that the x and y axes are 
no longer conjugate. On the other hand, the minimum can always be 
readied by successive searches along two conjugate directions. For 
example, consider 


°** [o]- • 

These vectors are conjugate since 


^ A u 2 = [ 1 0 ] 


5 I 



2 2 


5 

-3/2 5 


1 

2 

m « 



= 0 


and, as can be seen in Figure 6. 3, the minimum is reached by 

.2 

successive searches along u^ - [ 1 , 0 ] and u *1 5 , l] . 

Since the principal axes of a. quadratic function are orthogonal 
and also A -conjugate, one can always find the minimum by searching 
along the principal axes once only. The algorithm modifi cation by 
Brent is essentially based on this concept, i. e. , it is to find the 
principal axes of the function f ( or its quadratic approximation) and 
to search along the principal axes to obtain the minimum. 


6. 3 PRAXIS - A Modification of Powell' s Method 

Because of the deficiencies of Powell' 8 method, e.g., 

25 

the loss of linear independence and conjugacy, Brent developed a 
modified version called PRAXIS. The main modifications are: 
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1) A restart device is included to reset the search directions to 

a set of orthogonal, A-conjugate vectors after every n or n+ 1 iterations 
to insure the linear independence of the new search directions. These 
conjugate vectors are computed on the assumption that f is quadratic 
or is the quadratic approximation of the function to be minimized. If 
f is quadratic or if the quadratic approximation is good, then t-he new 
search directions are conjugate with respect to a matrix which is close 
to the Hessian matrix of f at the minimum. This resetting method 
will prevent the scheme from searching for a minimum in a subspace. 

2) A random step is inserted to enable the scheme to search for 
another initial point in each iteration if the most recent linear search 
has failed to improve the current approximati >n to the minimum. With 
this step in the scheme, the trouble noted by Zangwill will be avoided. 

For example, in Figure 6.1 if point B is chosen as the initial 
point, then Powell' s basic procedure will find C as the minimum and 
8 top, as noted by Zangwill. Powell' s modified procedure will retain 
the old search vectors as the new search direction for the next iteration; 
hence one more iteration is needed to reach the minimum. 

With the random step, the algorithm will replace point B by an 
arbitrary point in the space, say A' in Figure 6.1, after having failed 
to ot \in an improvement in the direction of u^ = (1, 0). This rules out 
the possibility of linearly dependent search directions. 

3) Discarding criterion. Powell' s modification proposed that 

the search direction should be discarded and replaced by one which 
maximizes ( det (V^ ) |, where 



Ml 


V. = ( u . A u.) 2 u.. 1< i ! n 
11 1 i — - 


(6.7) 


A: nxn matrix related to the quadratic approximation 
This discarding method may lead to the elimination of one of the mutually 
conjugate directions, in which case finite convergence for a quadratic 
function can no longer be assured. 

In PRAXIS the criterion described below is employed. It is 
essentially Powell 1 s criterion except an additional restriction is imposed 
to insure the finite convergence for a quadratic function property. 

The discarding criterion fcr PRAXIS is as follows: 

At K** 1 iteration with search direction u^ u^ and | det ( V |. • •. 

(a) For i = 1, . . . , n-k+1, take u. out of u^ t .... u^j^j » ‘‘Od compute 

T *1 

V T = f x " x ) A(x -x )] 2 (x -x ) 

J no no no 

(b) Compute |det(V .... V.^. V. +1> . . . . V n y Vj) | = D., i =1,..., n-k+1. 

(c) If no D. is larger than the value idet (V ,, ... V ) I , then no 

i 1 n o 

replacement for the search direction occurs. Otherwise go 
to (d). 

(d) For D = Max (D.) and D >1 det (V,, .... V ) | , renumber 

m i m 1 n o 

"l = "l \n-l' u m-l" u m = “m+1 Vl* V and 

U = X - X 

n n o. 

Thus, at the iteration, only one of u^, . . • » u n ^+1 permitted 
to be discarded. 

4) The linear search in PRAXIS is similar to Powell s procedure. 


It reduces the number of function evaluations considerably. For example, 



consider the linear search in the direction u, i. e. , minimize 

♦ (X) = f(x +A«), (6.8) 

o 

At the first iteration three function evaluations are reeded for a 

2 

quadratic cur*e - fit, say p (X) = aX + b A + c. The seond der- 

ivativeof P (X), i. e. , a. Is saved because it can be used in the next 

20 

iteration when this search direction is utilized again. Then, the 
approximation for the second derivative of p (X) is always available 
if a linear search in the direction u has already been performed or 
if u resulted from a singular value decomposition, which is the step 
to find the principal axis vectors in PRAXIS. Thus, only two additional 
function values are needed for the three constants a, b, c, where a = «j» "(0) 
after »he first iteration. 

6.4 Examples 

Zermelo' s problem is used to demonstrate the efficiency and 

17 

reliability of the algorithm. Long 1 s method is used to transform 

this variable time problem into a fixed final time problem. The equations 

of motion are 

x = V cos 0 

y = V 8 in 0 
♦ 

0 = u, |u | < 0. 5 (6.9) 

and the performance index is 

J =Ct f 2 tPj (x (t f ) - x f ) 2 +P 2 ( y (t £ ) - y f ) 2 +P 3 (0 (t £ ) - 0 f ) 2 (6.10) 

Performing Long 1 s transformation, i.e. , t = as, s« [0.1] . 



dt = ads 


f f :a 4,5:11 = 


Thus, 

x' = a V cos 0 
y* = a V sin 0 
0 1 = au 

a’ =0 ( 6 . 11 ) 

and 

J =Ca * fPj [ x (1) - x f ] 2 f P 2 [ y(l) = y f ] 2 +P 3 f flfl) - 0 f ] . 2 (6.12) 

To employ PRAXIS, the control u(s) must be discretized. Let u (s) = 

«J. * « [ 0.*j) ; U (s) = u 2> 8 « [ » 1 ,« 2 ] ; . • • ; u(s) = u q , s«[ s q ^ 1] . 

and consider the cos* J to be minimized as a function of the a + 1 

variables a, u_, . . , u , i.e., 

1 n 

J = Jf a,U|, a , . . , uj (6.13) 

The problem was then attacked with PRAXIS for three different initial 
estimates for a, e. , The results are summarized in Table 6.1. 



Table 6.1. Parameters and Results with PRAXIS. 


Constants: C = 1 

p l = p 2 = 10000, p 3 = 0 
n = 10, 

a™ = 0.2, i = 1, .... 10 

kj £-5 

* f = 4. 0, y £ = 5. 0 

imal Control: u. = . 5 for a « [ 0 , . 3] 



c. = 0. 0 for s { [ 0, \l] 






.734 

b 

.71 

G 


6.730 

6.7408 


6.71 3.999 
6.71 3.999 


4. 999 
4.$99 



The cost reduction versus number of linear searches for cases 


1 and 3 are shown in Figure 6. 5, while the various control profiles 

are shown in Figure 6.4. In addition, the cost for = 4 0 i® 

shown in Figure 6. 5. This trial converged to a local minimum and 

no longer improved the cost. Note that the trials w ith t f (0) = 1 and 

( 0 ) 

t £ = 10 converged to the neighborhood of the minimum cost rapidly. 
However, the control profiles, in a sense, oscillate about the optimal 


control. Of course, this behavior could be improved by assuming a 


more representative parameterized contro 





















= 1 sec. 



final control 
exac* solution 


(b) 4 sec. 
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(c) - 10 sec. 
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.8 1.9 |>.0 
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Figure 6.4. Control profiles using PRAXIS with various 
initial estimates of t^. 
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CHAPTER 7 


CONCLUSIONS AND RECOMMENDATIONS 

7.1 Summary 

Computer programs ter shuttle trajectory optimization have been 
developed and delivered to NASA- JSC. Cne of the programs contains 
the function -space gradient, conjugate -gradient, Davidon, and Broyden 
algorithms for the ascent problem. The PILAXIS parameter optimization 
scheme has been integrated into the NASA-JSC PEACE parameter 
optimisation program. 

7.2 Conclusions and Recommendations 

1. ) The function-space Broyden and Davidon methods performed 
appreciably better on the shuttle ascent problem than the gradient and 
conjugate -gradient algorithms, with Broyden slightly better than 
Davidon. Both control and state variable inequality constraints were 
included in the formulation with the control constraints handled directly 
while the state variable constraints were included with penalty functions. 

2. ) The storage problems associated with function -space Da> r don- 
type techniques have br n overcome. Although considerable storage 

is necessary for the computation of inner products, the storage need 
not be in fast memory. On the University oi Michigan computer the 
storage problem is handled very easily by a disk file storage system. 

The programs for use on the NASA-JSC computer require modifications 
for drum storage. 
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3. ) The University of Michigan computer is a time -sharing system 
with an interactive -graphics capability. This capability accelerates 
considerably the time required to converge a large-scale optimization 
problem. For exnmple, standard operation with the NASA-JSC PEACE 
parameter optimization program (when a large number of parameters 

is involved) usually involves: making a single computer run daily, 
analysis of the result, adjustment of parameters (usually penalty 
coefficients), and resubmission of the program. This means that the 
analyst must stop-and-start on the same problem many times, and the 
process is a somewhat inefficient use of the analyst 1 s time. With a 
tirue -shared, interactive graphics capability, the analyst can stay with 
the problem continuously for longer periods of time with the result 
being: less total computer time, less total human effort, more physical 
knowledge of the problem, and more rapid solution of the problem. 

Thus, it is recommended that MPAD consider the use of interactive- 
graphics terminals in the solution of large-scale trajectory optimization 
and mission analysis problems. 

4. ) Previous investigators have noted difficulties in solving 

variable final-time trajectory optimization problems with accelerated- 

gradient methods. In Chapter 5 heuristic ar^-ments and simulations 

(0) * 

indicate that the initial estimate of t^ is critical, and t^ < t^ improves 
the convergence rate considerably. 
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5. ) Due to budget limitations, the PRAXIS algorithm could 
not be simulated on realistic shuttle trajectory optimisation problems . 
The worth of this algorithm will be determined by NASA- JSC personnel. 
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Appendix A 


Dynamical Equations of Motion: First Stage 
As is customary in trajectory optimization the vehicle is modeled 
as a point mass. It is furth* r assumed that the thrust, aerodynamic, 
and gravitational forces act through the center of mass. By Newton 1 s 
Second Law. 

Zv F = m r , (A. 1) 

* * 

where r is measured in an inertial coordinate system. Since numerical 
integration is desired in the first stage, consider Figure A. 1. The 
acceleration of r is: 


Z 



r = R +« xp +wx(u xp)+P ROT Uwxp rqT ( a -2) 


Consider two coordinate systems fixed at the center of the earth, 
one of which rotates with the earth and the other inertial. Then, 

i) R = 0, since both coordinate systems are fixed at the same point. 

ii) 55 - constant <3 = 0, since rototation of earth about its axis 

is constant. 

iii) r = p; follows from r = R + p and i) R = 0. 
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Thus the acceleration is 
■ * •# 

r = r x(w x r) x - f \ . 3 ) 

rot Ro i 

Now consider the "wo sphe~*.ca' coordinate ^v>ter»*s centered at the 
center of the earth shown Ln Figure A. 2- By definition _ 



Non-Rotating rotating about z axis 

Figure A. 2. First-Stage Coordinate System. 

Since £2 is along the z-axis, 

£2 = £2 [Cos 9 e - Sin 9 e } , 


which implie^.__ 


£2 x r = 


r P 

ft Cos P~ ft Sir 

r n 


= rft Fin P 


0 


SI x (£2 x r) = 


SI Cos 9 -ft Sin 0 


0 


0 


= - [ r ft ^ Sin^ Q ] e 


0 

r f2 Sin 0 . 

O 

r ft^ Sin 0 Cos 9] e 
J 0 


r 

" I 
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2 fi x r 


ROT 


e 


2 A Cos 0 


V - V 

0 ♦ 

- 2 Cl Sin 0 0 

r 0 r$Sin 0 


= - [ 2 r £2 <|» Sin 0] e - f 2 r Cl 4> Sin 


+ [ 2 r Cl 0 Cos 0 + 2 r Cl Sin 0] 

yet T = Thrust Force = T e + T e + T e 

r r to 0 <p 9 

A = Aerodynamic Force = A^e^ + A^ e^ + 

— — yn k — 

G - Gravitational Force = r e 


■<!> 


Then, upon substitution into Eq. (A. 1) 

• 2 *2 2 
r - r 0 - Sin 0 

-rs/ Sin^ 0 - 2 r £2<)> Sin ^ 0 

1 T _ . k , 

=- [ T + A - m — ] 
m 1 r r 2 J 

r 

• • * • ^ 
r0+2 r0-r<|> Sin 0 Cos 0 

• rl/ Sin 0 Cos 0 -2 r ft <|> Sin 0 Cos 0 

= - [ T + A ft ] 
m 0 0 1 

m « • * • 

r <J> Sin 0 + Z r <l> Sin 0 + 2 r 0 4* Cos 0 * 

* 

12 rf!0 Cos 0 f 2 r Cl Sin 0 


= - [ T , + A . ] 
m 1 q> <}> J 


Define: 


r = v 


0 = v 


0 


<j> = v 


<t> 

r Sin 0 


0 Cos 0 ] 


(A. 4) 


(A. 5) 


(A. 6) 
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m = mass flow rate 


F = A + T 
r r r 


F = A + T 

e e e 


F 4> = A 4> +T «t» 


Since is an ignorable coordinate, the 4> - equation is neglected. Then, 
with the following state-control definitions 

*1 = r - R o' X z = 9 ■ X 3 = V x 4 = v $ ' x 5 = V *6 =m - “ * ,m-1 


the equations of motion are: 

(r) \ = X 3 

<e » Mv* 7 *,) 

X + X 2 

(u , i = i +(x + R Jo 2 Sin‘ * 

3 + R o ) <V R 0 I 1 J 2 


+ 2 fi x Sin 
5 


,V) k 4 * S " X 2 CoS X 2 


+ 2 ft x_ Cos x_ + — 

5 2 x 6 


X 3 X 5 

(w) *5 = -(V^T 


X 4 X 5 

Cj + lytan’x^ 2 nx 4 Cos x 2 


(A. 8) 


- 2 x 0 ft Sin + — 
3 2 x 6 


(mass ) = - u 
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where 


*1 

X Z 


- altitude above earth 

= e 


x, = v -velocity in e direction 
3 r 7 r 


x = v - velocity in e direction 
4 0 0 

Xg = v^ - velocity in e^ direction 


x^ = mass of vehicle 
u - [m | - mass flow rate 
R = radius of earth 


o 

£2 - angular velocity of earth 
k = gravitational constant of earth 



Appendix B 


Atmosphere and C 


D 


Models 


The 1963 Patrick Atmosphere model was used. Pressure and density 
ratios, and speed of sound data were obtained from Ref. 26 and curve 
fitted as functions of altitude according to the equations 

P ^SL~ “P < a 0 + a i x + ■ • 1 • + a n x J 

pfPgjj e*P ( b 0 + b j * + • • • • + b i3 x 13 ) 

13 

a = exp (c+c, x+....+cx ) 

o 1 13 

altitude (ft) - 200, 000 

x = 100, 000 

The coefficients a., b^, c. are given in Ref. 27 

The data in Table B.l were used to define the drag model. The drag 
force is given by 



V*C d A 


For Mach numbers between those in the data table, interpolation by 
piecewise cubic splines was used. 
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Mach No. (M) 

C D 


.028 

0.2 

.028 

0.4 

.028 

0.5 

.029 

0.6 

.032 

0.8 

.058 

1.0 

.110 

1.1 

.121 

1.2 

.123 

1.5 

.121 

1.75 

. 115 

2.0 

.106 

2. 4 

.088 

3.0 

. 066 

3.5 

.055 

4.0 

.047 

5.0 

.039 

6.0 

.031 

7.0 

.025 

8.0 

.021 

9.0 

.019 

0.0 

.019 


Table B.l Drag Coefficient 








Appendix C 


Transformation Equations 

At the end of first -stage burn we wish to determine both the 
inclination of the orbital plane and the initial conditions for the second- 
stage, non-rotating polar coordinate system. The plane of the orbit is 
determined at first-stage burnout because there is no out-of-plane 
thrusting in the second-stage. Consider the conditions at fir?t-sta 
burnout: 

r = r e 

r 


v s + fix r = v e + v e + v e + r OSin9 e 

ROT r r 0 9 9 9 9 


= v e + v n e„ +( v, + r fi Sin 9) e , 
r r 8 0$ 4> 

Consider the new polar system ( 0~, v”\ v ) : 

r B 



Figure C. 1. Second-Stage Polar Coordinate Sy stem, 
since Ivi = f 7 " ' ,Z 


(C.l) 


ince |v| = /v “ + v«“ + ( v +rfi Sin 0 )“ and the radial components of 
'It “ <j> 

the velu^Uy are the same in each system, the velocity transformation is 


a/ 

V -V 

r r 


'v_ 


=Jv 2 + ( 

N 9 


(v + r fi Sin 9 ) 


Thus, in state notation 


x i = x l 


Y r v 

2 3 




= Jx 4 Z + + C ^ + R q ) fi Sin x 2 ] T 


(C.2) 


(C.3) 
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where x,= r-R . xlw" , x*„ = 91 . Note that the mass will change by the 
1 o 2 r» 3 0 

amount of structure discarded. 

To obtain the inclination, consider the following unit ver r which 

is perpendicular to the plane of the orbit: -> 

_ _ _ C r C 0 % 

N = r x V = r 0 0 

v v (v . + r SJ Sin 0 ) 
r 0 <|> 

= rV 0 % ‘ r(v * +rJlSinG ) e 0 


S N = 1W] 


7^e +( % + rQSin9) 

(v + r £1 Sin 0 ) 

♦ 

/ yZ 0 + (v 4> + rnsin0)Z 


Let k ~ unit vector along the axis of rotation of the earth. Then, 

k ~ Cos 0 e - Sin 0 e. 

r 0 

The relation between e. T , k , and the inclination, ft, is shown in Figure C.2, 

N 


c n r 


ft = inclination angle 


Figure C.2. Orbit Inclination Geometry. 


Thus, 


v , + r £1 Sin 0 ) Sin 0 


Cos i = e • k = 
N 


V- + ( v . + r £1 Sin 0 ) 
0 9 



o*- in ^trite notation 


' x + (x.+R ) ft Sin x 1 Sin x 

Dio 2' 2 i 

(C.4i 

f* y 2 " — 

!>: +■ f x_ + (x„ + R ) ft Sin x 1 |2 

14 ! d 1 2 J 



Denote the transformation equations (C* 3' by 

^ = <i ( x). (C.5) 

The partial derivatives of g with respect to the x. are defined by 

i) v*i = 




i = 2,6 



Let Xj.-(Xj+R o )ft Sin = ( ) 

! * 4 2 • ' > 2 1 = f 1 

Thus 



ag 3 2 - 
Sx, = 2 1 


-1 

] 2 2 ( ) ft Sin x^ 

-1 

] 2 2 ( )(x. +R ) Cos x 

i o c 
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*g 


3 



= « 






3 


Bx 


5 



li 


2 ( 


) 



Appendix D 

User s Guide for PRAXIS 

PRAXIS determines the local minimum of a scalar function 
which need not be differentiable. Double precision is necessary for 
all floating point variables. An EXTERNAL statement for ►he function 
to be minimized is needed in the program which calls PRAXIS. 

However, the gradient of the function is not required. 

Usage of PRAXIS 

CALL PRAXIS (TO, HO, N, IPRIN, X, F, FMIN) 

Description of parameters: 

F : Function to be minimized 

MACHUP: A machine precision parameter furnished in the 

program; it is about 2.22 x 10 ^ on the IBM 369. 

TO: A tolerence for the stopping criterion; the program 

stops searching for the minimum if 

| | x l -x l + 1 j |£ MACHUP | | x I+1 1 | + TO 

HO: Maximum step-size. To assure fast convergence, 

HO should be about the maximum distance from vhe 
initial guess to the minimum. 

N; The number of dependent variables, i. e. , the dimension of 

x (N should not be less than two). 

IPRIN: An integer for controlling the printing of numerical 

results. 


12 3 



1?.4 


IPRIN 

IPRIN 


IPRIN 


IPRIN 


IPR7N 


F(X, N) 
FMIN : 

Out put v'o riable s . 

JLiVlIN: 

E VAL3: 
MIN F: 


0. Nothing is printed by PRAXIS. 

1. Value of F is printed after every X + I or X 2 
linear minimizations. Final x is printed. If 
X 1, intermediate x is printed also 

2. The scale factors and the principal values of the 
approximating quadratic form are also printed. 

3. The values of x after every few linear minimizations 
are printed also. 

4. All available and relevent values are printed. 

An N dimensional vector. Initial guess of minimum 

is placed here to start the program. Final estimate 

of X is returned to here. 

A REAL * 8 function to minimized. A declared 

EXTERNAL is necessary in the calling program. 

The final value of F obtained. 

Number of linear minimizations. 

Number of function evaluations- 

Function value at LMIN th linear minization 


Example of use 

IMPLICIT REAL * 8 (A-H, ff -Z) 
DIMENSION X (2) 

EXTERNAL BANANA 


ro = 1. D-5 
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N = 2 

X(l) =-1.2 DO 
X (2) = 1. DO 
HO = 2. 0 
IPRIN = 1 

CALL PRAXIS (TO, HO, N, IPRIN, X, BANANA, FMIN) 

PRINT FMIN 

1 FORMAT ( ' FIv£IN =' , D 25.15) 

END 

C. . . Fivuction to be minimized 

REAL FUNCTION BANANA (X,N) 

IMPLICIT REAL * 8 (A-H, $ -Z) 

DIMENSION X (N) 

-A- 

BANANA = 100. DO * (X (2) - X (1) ** 2} 2 + (l.DO- X(l* ) ** 2 

C . . .NOTE, THERE ARE NO DERIVATIVES OF BANANA 

RETURN 


END 









